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FOREWORD 



This final report describing the formulation of the Program 
o Optimize Simulated Trajectories (POST) is provided in accord- 
ince with Part IV of NASA Contract NAS1-13611. The report is 
.resented in three volumes as follows: 

Volume I - Program to Optimize Simulated Trajectories - 
Formulation Manual; 

Volum II - Program to Optimize Simulated Trajectories 
Utilization Manual; 


Volume III - Program to Optimize Simulated Trajectories 
Programmer's Manual. 

This work was conducted under the direction of Joseph Rehder of 
the Space Systems Division, National Aeronautics and Space 
Administration, Langley Research Center. 



11 


^ M r iwmjm WflWW .W.W , AVI 1 


/ 


CONTENTS 


SUMMARY 


INTRODUCTION 


LIST OF SYMBOLS AND ABBREVIATIONS 


COORDINATE SYSTEMS 

Coordinate System Definitions 

Attitude Angles 

Transformations 


PLANET MODEL , . . 
Oblate Spheroid . • 
Gravitational Model 
Atmosphere Models . 


VEHICLE MODEL 

Mass Properties Model . 
Propulsion Calculations 
Aerodynamic Calculations 
Aeroheating Calculations 
Steering Model • . • . 


TRAJECTORY SIMULATION . 

Events/Phases 

Translational Equations 

Static Trim 

Integration Methods . • 
Launch Options .... 


AUXILIARY CALCULATIONS 
Conic Calculations . , 
Rangft Calculations • , 






} 


I 


* 4 , 



i 

l 


Auxiliary Position and Velocity Calculations 
Auxiliary Attitude Calculations «••••• 

Tracking Data 

Analytic Impact Calculations 

Velocity Losses 

Velocity Margins 

Sun-Shadow Calculations 

Multiple Vehicle 


VI1-7 
VII-10 
VII-12 
VII-14 
VH-16 
VII-16 
VIl-17 
VII-18 
and 
VII -19 


VIII. TARGETING AND OPTIMIZATION 

Problem Formulation 

Fundamental Concepts and Nomenclature 

Direction o£ Search 

Step-Size Calculation 

One-Dimensional Minimization .... 
Algorithm Macrologic 


VII1-1 
VILI-1 
VII 1-4 
VIIL-10 
VIII-16 
VIII-19 
VIII-24 
thru 
VIII-27 


...i 


i 


t 


IX. REFERENCES 

APPENDIX A 

DECOMPOSITION BY PARTITIONING INTO FULL-RANK 
SUBPROBLEMS 


IX-1 


A-l 
thru 
A- 7 


Figures 


1-1 Program Macrologic 

III-l Coordinate Systems 

111-2 Launch Frame 

II1-3 Body Frames 

LII-4 Orbital Parameters ....... 

III-5 Inertial Euler Angles 

III-6 Relative Euler Angles 

III- 7 Aerodynamic Angles . • « • - • • 

LII-8 Inertial Aerodynamic Angles . • 

TII-9 Matrix Transformations 

IV- 1 Oblate Planet 

V- l Engine Gimbal Angles 

V-2 Aerodynamic Angles 

V- 3 Body Rates * * 

VI- 1 Event Sequence Setup 

Vi-2 Illustration of Time-to-Go Logic 

VI-3 Moments in Pitch 

VI-4 Moments in Yev 


1-4 
I1I-2 
111-2 
I1I-3 
1 11-4 
111-5 
111-5 
HI-6 
I1I-6 
m-7 

IV- 3 
V-5 
V-6 

V- 15 

VI- 1 
Vl-4 
Vl-9 

VI-10 


iv 


Iff" 


Down range and Crossrarge Angles VII-4 

Radar Tracking Schematic VI I -12 

Relative Geometry between Active and 

Target Vehicles VII-19 

Direction of Search in the Independent- 

Variable Space VIII-10 

Illustration of Minimum-Norm Constraint 

Direction for n * 2 < m * 3 VI1I-11 

_ a 

Illustration of Newton-Raphson Constraint, 

Direction n « m • 3 VIII-12 

a 

Illustration of Least-Squares Constraints 

Direction for n =4 > m “ 3 VIII-13 

a 

Direction of Negative-Projected Gradient 

for n * 1 and m = 3 VIII-15 

a 

Properties of Estimated Net Cost Function VIII-18 

Macroiogic of Projected Gradient Algorithm VIII-25 

Complete PGA Iteration, Construing of Optimisation 
Step Followed by Constraint Step for n^ ■ 1 

and m - 3 VIII-27 

Tables 


Typical Applications of POST 

1962 I). 5. Standard Atmosphere Profile • . . < 
Derived Coefficients for the 1963 Patrick AFB 
Atmosphere Model 


1963 Patrick AFB Molecular Temperature Profile 

and Gradient Profile 

Runge-Kutta Coefficients (4th Order) 

Runge-Kutta Coefficients (8th Order) 

Laplace f and g Series 




JNW F - w f" 

> r i 



e 

r 

r 

r 

k- 

r • \ 


*■ / 

t * 1 



i 

| 



r 1 


i 

i 


FINAL REPORT 

PROGRAM TO OPTIMIZE SIMULATED TRAJECTORIES (POST) 

VOLUME I - FORMULATION MANUAL 

By G. L. Brauer, D. E. Cornick, A. R. Habeger, 

F. M. Petersen, and R. Stevenson 
Martin Marietta Corporation 


SUMMARY 


This report documents the equations and the numerical tech- 
niques used in the Program to Optimize Simulated Trajectories 
(POST). 

POST, a generalized point mass, discrete parameter targeting 
and optimization program, provides the capability to target and 
optimize point mass trajectories for a powered or unpowered 
vehicle operating near a rotating oblate planet. POST has been 
used successfully to solve a wide variety of atmospheric flight 
mechanics and orbital transfer problems. The generality of The 
program is evidenced by its f!-pha$e simulation capability, which 
features generalized planet and vehicle models. This flexible 
simulation capability is augmented by an efficient discrete 
parameter optimization capability that includes equality and 
inequality constraints. 

POST was originally written in FORTRAN IV for the CDC 6000 
series computers. However, it is also operational on the IBM 
370 and the UNI VAC 1103 computers. 

Other volumes In the final report are: 

Volume II - Utilization Manual - Documents Information 
pertinent to users of the program. It describes the 
input required and output available for each of the 
trajectory and targeting/optimization options. 

Volume III - Programers Manual - Documents the program 
structure and logic, subroutine descriptions, and other 
pertinent programing Information. 
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I. INTRODUCTION 


POST is a general purpose FORTRAN program for simulating and 
optimizing point mass trajectories of aerospace type vehicles. 

The program can be used to solve a wide variety of performance 
and mission analysis problems for atmospheric and orbital vehicles. 
For example, typical applications of POST are outlined in Table 
1 - 1 . 


Table 1-1, - Typical Applications of POST 


Typ« of 
Him* lot) 


! 

Optlalxatlon 

Variables 

Typ veal Constraint a j 

CPU Tine Required 

Typo of Vehicle 

2qu* llty 

Inequality 

nln 

Ascent to 

N«ir*E«rth 

Orbit 

Titan IIIC A ME, Space 
Shuttle, Single State 
to Orbit (VTO end HTO) 

Payload, Weight at 
Burnout Fuel, Burn tins. 
Ideal Velocity, 

Initial Weight 

Radius 

Plight Path Angle 
Velocity 

Dynanlc Pressure 
Accelerations 

2 - 20 

Ascent to 
Synchronous 
Equatorial 
Orbit 

Titea IIIC, Shuttle/ 
Tug 

Payload 

Apoge* 

Perigee 

Inclination 

Angle of Attach 
Pitch Rates 

3-50 

As cone 
Abort 

Spec# Shuttle 

Abort Interval 

Landing Site 
Longitude and 
Latitude 

Acceleration 
Dynanlc Praaaurn 

2-5 

ICBM 

ballistic 

His* Urn 

Titea 11. 
Want—— I A U» 
Ufa— ard 

Payload 
Mac Distance 

Latitude 
Longitude 
Cross rang* 
Ho— range 

Flight Path 
Angle at Bntry 
Acceleration 
during Entry 

2 - 2Q 

Koontry 

Space Shuttle i X-24C, 
Slogl* ftAf« to Orbit 

Heat Rata 

Total Heat 

Crosarang* 

1 

latitude 
| Longitude 
i Cro* avenge 
Do— range 

Hast Rate 

Acceleration 

3-15 

Orbital 

Maneuvers 

Trans tags. 

Space Tug, 1US, 
Solar Electrical 
Propulsion 

Payload 

P—1 

! 

fcadlue 
Ve* K lty 

Path Angle 

Ar^—an t of 
Perigee Period 
/ipog— , Perigee 

Attitude Angles 
Perigee Altitude 

0.5 - 10 

Aircraft 

Nrfomaci 

X-241 aad C, Space 
Shuttle SObecale, 
Subsonic 

Jet Cruiae, Hypersonic 

Boabara and • 

Inter cap tore 

Mach Mu— er 
Cruiae Tine 
Payload 

Do— range 
Croaa range 
Dynsnlc Pressure 
Velocity and 
Mach Altitude 

Dynanlc Pressure 

Dynanlc Pressure 
at Mat Altitude 

0.1 - 5 


One of the key features of POST is an easy to use NAMELIST- 
type input procedure. This feature significantly reduces input 
deck set-up time (and costs) for studies that require the normal 
large amount of input data. In addition, the general applicability 
of POST is further enhanced by a general-purpose discrete parameter 
targeting and optimisation capability. This capability can be 
used to solve a broad spectrum of problems related to the per- 
formance characteristics of aerospace vehicles. 

The basic simulation flexibility is achieved by decomposing 
the trajectory into a logical sequence of simulation segments. 
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These trajectory segments , referred to as phases, enable the tra- 
jectory analyst to model both the physical and the nonphysical 
aspects of the simulation accurately and efficiently. By segment- 
ing the mission into phases, each phase can be modeled and simu- 
lated in a manner most appropriate to that particular flight regime. 

For example, the planet model, the vehicle model, and the simula- 
tion options can be changed in any phase to be compatible with the 
level of detail required In that phase. 

t 

Every computational routine in the program can be categorized *• t 

according to five basic functional elements. These elements are: 
the planet model, the vehicle model, the trajectory simulation 
model, the auxiliary calculations module, and the targeting and 
optimization module. The planet model is composed of an oblate 
spheroid model, a gravitational model, an atmosphere model, and 
a winds model. These models define the environment in diich the 
vehicle operates. The vehicle model comprises mass properties, 
propulsion, aerodynamics and aeroheating and a navigation and 
guidance model. These models define the basic vehicle simulation 
characteristics. The trajectory simulation models are the event- 
sequencing module that controls the program cycling, table inter- 
polation routines, and several standard numerical integration 
techniques. These models are used in numerically solving the 
translational and rotational equations of motion. The auxiliary 
calculations module provides for a wide variety of output calcu- 
lations. For example, conic parameters, range calculations, and 
tracking data are among the many output variables computed. The 
targeting and optimization module provides a general discrete 
parameter iteration capability. The user can select the ootimiza- 
tion variable, the dependent variables, and the independent vari- 
ables from a list of more than 400 program variables. An accel- 
erated projected gradient algorithm is used as the basic optimiza- 
tion technique. This algorithm is a combination of Rosen's pro- 
jection method for nonlinear programming and Davidon's variable 
metric method for un cons train ted optimization. In the targeting 
mode, the minimum norm algorithm is used to satisfy the trajectory 
constraints. The cost and constraint gradients required by these 
algorithms are computed as first differences calculated from 
perturbed trajectories. To reduce the costs of calculating 
numerical sensitivities, only that portion of the trajectory in- 
fluenced by any particular independent variable is reintegrated 
on the perturbed runs. This feature saves a significant amount 
of computer time when targeting and optimization is performed. 

POST is operational on several computer systems as described 
in the tabulation. 
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Location 

Computer 

Operating System 

Martin Marietta Corporation 
Denver, Colorado 

CDC 6400, 6500 

SCOPE 3.4.1 

Martin Marietta Corporation 
Michoud, Louisiana 

UN IV AC 1110 

EXEC 8 

Langley Research Center 
Hampton, Virginia 

CDC 6600 

SCOPE 3.2 

Johnson Spacecraft Center 
Houston, Texas 

UN I VAC 1108 

EXEC 8 

Goddard Spaceflight Center 
Greenbelt, Maryland 

IBM 370-192 

OS 

Marshall Spaceflight Center 
Huntsville, Alabama 

UNI VAC 1108 

EXEC 2 


Basic program macrologic is outlined in figure I— 1, which 
illustrates the linkage between the simulation and the Iteration 
modules . 
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LIST OF SYMBOLS AND ABBREVIATIONS 



Math symbol 

Internal Fortran symbol 

Definition 

a 

SEMJAX 

semimajor axis, m (ft) 

A,A 

— 

comr. orient of radius vector 
perpendicular to Sun vector , 
m (ft) 

— AB * ( A AXB» A AYB y A AZb) 

— 

aerodynamic acceleration 
in the body frame, mps 2 
(fps 2 ) 

[AB] 

AB(I) 

matrix transformation from 
the A-frame to the B-frame 

A 

c 

— 

centrifugal acceleration, 
mps 2 ( f ps 2 ) 


AR 

nozzle exit area of each 
rocket engine, m 2 (ft 2 ) 


AHORIZ 

horizontal acceleration , 
mps 2 (fps 2 ) 

kt 

AXIT, AYIT, AZIT 

acderation of target 
vehicle in the (ECI) frame, 
mps 2 (fps 2 ) 

A i’ *i 

a m* A MP» ^ 

AMXB, AMYB, AMZB 

constants 

totcl aerodynamic moment 
about the roll, pitch, yaw 
exes, l»-m (ft-lb) 

N 

A(I) 

Davidon deflection matrix 
component 

A s 

ASM 

total sensed acceleration, 
mps 2 (fps 2 ) 

^SB " (^XB* A SYB* A SZb) 

AXB, AYB, AZB 

total sensed acceleration 
in the body frame, spa 2 

( tv* 2 ) 

^SI “ (^XI* ^YI* ^Zl) 

AS XI, ASYI , ASZI 

total sensed acceleration 
in the Inertial frame, 
■ps 2 (fps 2 ) 

^TB " ( A TXB* ^TYB* K . IU ) 

— 

thrust acceleration in the 
body frame, ape 2 (fps 2 ) 





^TB " (SxB* ^TYB* Sa) 
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Math symbol 

Internal Fortran symbol 

Def ini tlon 

s 

AVERT 

vertical acceleration, 
mps 2 (fpa 2 ) 


AZL 

azimuth of the z axis, 
rad (deg) 

^.1’ A ZR* ^ 

AZVELI, AZVELR, 
AZVELA 

azimuth of the inertial, 
relative, and atmospheric 
relative velocity vectors, 
rad (deg) 


AZWT 

wind azimuth, rad (deg) 

a zt 

TKAZMI 

azimuth of the slant range 
vector to the tracking 
station, rad (deg) 

B i 

— 

boundary for i^ con- 
straint 

B 

tl 

B(l) 

Davidon deflection matrix 

B(R) 

— 

boundary of region R 

B(u) 

— 

local boundary hypet 
surface 

C A’ S* °N 

CA, CY, CN 

axial, side force, ami 
normal aerodynamic force 
coefficients 

A o o o 

— 

component of C^, Cy, Cj, 

that is not multiplied 
by a mnemonic multiplier 

V °L 

CD, CL 

drag and lift coefficients 

V \> 

— 

drag and lift coefficient 
components that are not 
multiplied by a mnemonic 
variable 
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Math s ymbol 


Definition 


C(u) 

D 

E 

IE] 



CM, CW 


CS 

E(I) 

DRAG 

ECCAN 

ECCEN 

E0(I) 

E(I) 


WE(1) 

F 0PTVAR 

f 


^AB “ 1 

|F 

( AXB* 

f ayb* 

f azb 

| FAXB, 

FAYB, 

FAZfe 

II 

If 

^ TXB* 

f tyb* 

f tzb] 

| FTXB, 

FTYB, 

FTZB 


[GA] GA(I) 

GHA, GHAS GHA, GHAS 


pitch and yaw moment co- 
efficients 

speed of sound, mps (fps) 

constraint functions 
aerodynamic drag, H (lb) 
eccentric anomaly 
Euler parameter matrix 
eccentricity 
Euler parameters 

active constraint error 
vector. 

weighted error vector 
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optimization function $ 

s 

nonlinear vector-valued 1 

function ! 

I 

aerodynamic forces in the < 

body frame, M (lb) 

thrujt forces in the body 
frame, N (lb) 

matrix transformation from 
the G-frame to the A-frame 

Greenwich hour angle and 
Greenwich hour angle of Sun, 
rad (deg) 


S " ( G xi* Si* G zi) 


GXI, GY1, GZI total gravitational 

acceleration in the ECI-- 
frame, mps 2 (fpa 2 ) 
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Math symbol 


Internal Fortran symbol Definition 


*n 


H 

h 


DG(I) difference In the gradient 

vector VF between the cur- 
rent and previous Itera- 
tion 

gravitational constant 

ALTIT0 oblate altitude, a (ft) 


- “ (\l* Nl* ^i) 


k 


h 

P 


H 


B 


h 


c 


ANGM0M 


angular momentum, mps 2 
(fps 2 ) 


ALTA, ALTP 
HB 


P2 


altitude of apogee and 
perigee, km (n mi) 

base altitude used in 
atmospheric calculations, 
m (ft) 

constraint function 


H 


g 


HT 


geopotential altitude, m 
(ft) 



heating ratios 


TRKHTI altitude of tracker, m 

(ft) 

P1NET estimated r.et cost func- 

tion 


i 

[IB] 


[IG] 


INC 

relative-frame orbital 
inclination, rad (deg) 

IB(I) 

matrix transformation from 
the ECI-fraiu to the body 
frame 

IG(I) 

matrix transformation from 
the ECI-frame to the geo- 
graphic frame 
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Math symbol 

Internal Fortran symbol 

Definition 

HU 

IL(I) 

matrix transformation from 
the ECI-frame to the 
launch frame 

[IP] 

IP(I) 

matrix transformation from 
the ECI-frame to the 
planet frame 

I 

sp 

I3PV 

rocket specific impulse, 

%J 

[J] 

AC0B(J) 

constraint Jacobian matrix 

l I ] 

2’ 3’ 4 

J2 , J3 , J4 

gravitational constants 

k - (k x , k 2 , k 3 , k 4 ) 

— 

Runge-Kutta constants 

K i 

— 

constants 

L 

LIFT 

aerodynamic lift, N (lb) 

[LB] 

LB(1) 

matrix transformation 
from the launch frame to 
the body frame 

i 

LREF 

aerodynamic reference 
length, m (ft) 

M 

MACH 

Mach number 

M 

MEAN 

mean anomaly, rad (dag) 

M 

— 

pitch and yaw moment equa- 
tions 

in 

MASS 

vehicle mass, kg (slug) 

M f 

— 

mnemonic table multiplier 
for table f 

n 

a 

NAC 

number of active con- 
straints 
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a 

c 

NDEPV 

number of constraints 

[pi, ih 

PR0J (I) 

projection operators used 
in the projected gradient 
method 

p(h) 

PRES 

atmospheric pressure, 
N/m 2 (psf) 

p i 

PI 

weighted optimization 
variable 

P 2 

P2 

weighted constraint error 
function 

Q 

TLHEAT 

total heat, J/m 2 (Btu/ft 2 ) 

q 

DYNP 

dynamic pressure, N/m 2 
(lb/ft 2 ) 

• * 
^lam’ ^turb 

HRATRT, HTURB 

laminar and turbulent heat 
rate, W/m 2 /s (Etu/ft 2 /s) 

Q(o) , q(u) 

— 

linear manifold and its 
orthogonal complement 

r a 

RTASC 

right ascension of out- 
going asymptote, rad (deg) 

r 

a 

AP0BAO 

apogee radius, m (ft) 

[R3] 

RB(I) 

matrix transformation from 
the body reference frame 
to the body frame 

*D 

DPRNG1 

dot-product range, km 
(n ml) 

R E’ *P 

RE, RP 

equatorial and polar 
radius, a (f.) 

R» 

RN 

nose radius, a (ft) 
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Math symbol 

Internal Fortran symbol 

Definition 

*W 

REYN0 

Reynolds number 

r, - ( V y,, tj) 

XI, YI, ZI 

inertial radius vector 
from center of planet to 
the vehicle, m (ft) 

r i 

GCRAD 

geocentric radius, m (ft) 

-It 

XIT, YIT, ZIT 

radius vector to target 
vehicle, m (ft) 

r 

P 

PGERAli 

perigee radius, m (ft) 

R 

s 

RS 

radius to oblate surface, 
m (ft) 

^SR 

— 

slant range vector, m (ft) 

^SRG 

— 

slant range vector in geo- 
graphic frame, m (ft) 

^TR 

— 

radius vector to tracking 
station, m (ft) 

^VE 

XVE, YVE, ZVE 

Radius vector of vehicle in 
vernal equinox system, m (ft) 

8_ 

S(I) 

direction of search 

C 

£ 

— 

direction of seerch to 
satisfy the constraints 

s 

SHADF 

shadow function, m (ft) 

^toss^ 

SL0SIJ 

space losses for tracking 
stations, dB 

S r.f 

SREF 

aerodynamic reference 
area, n 2 (ft 2 ) 

c° 

— 

direction of seerch for 
optimisation 

X 

ATEM 

atmospheric teaperature, 
•K CF) 

c 

TIME 

time, s 

V* 

— 

Jet engine thrust, N (lb) 


II-7 


TTJCB 


Math symbol 
T M1P 

T M2P 

t miy 

T M2Y 


T°(y) 


R 



vac 


SP 


T 


TP 


U 


u 


Internal 


ran svmool 


TTMZB 


TMYB 


TMZB 


THRUST 


TVAC 

TIMSP 

TIMTP 

U(I) 


Definition 

total thrust moment for 
nontrlmlng engines in the 
pitch plane, body axis 
system, N-m (ft-lb) 

total thrust moment for 
trimming engines in the 
pitch plane f body axis 
system, H-n (ft-lb) 

total thrust moment in the 
yaw plane for nontrinming 
engines, body axis system, 
N-m (ft-lb) 

total thrust moment in the 
yaw plane for the trimming 
engines, body axis system, 
N-m (ft-lb) 

denotes n^ 1 order table 
interpolation on the 
variable y 

total rocket thrust for 
all engines, N (lb) 

total resultant rocket 
thrust for engine 1, N 
(lb) 

vacuum thrust for rocket 
engines, N (lb) 

time since perigee pas- 
sage, s 

time to next perigee pas- 
sage, s 

gravitational potential 
function 

independent variable 
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Math symbol 


on 


^AB ( U B* V B* w b) UB, VB, WI components of the atmos- 

pheric relative velocity 
vector expressed in the 
body frame, mps (fps) 


-&I 



unit vector along the 
radius vector 

u 

—B 


X<= ,T E, YSUE, ZSUE 

unit vector in Sun direction 
in the vernal equinox system 

hi 


XSI, YSI , ZSI 

unit vector in Sun direction 
in the ECI system 

^1 



unit vector along the 
velocity vector 

All 


DU(I) 

change in the independent 
variables 

V 

a 


APVEL 

inertial velocity at 
apogee, mps (fps) 

-AG 


UA, VA, WA 

atmospheric relative 
velocity in the G-frame, 
mps (fps) 

hi 


VAXI, VAYI, VAZI 

atmospheric relative 
velocity vector in the 
inertial frame, mps (fps) 

h" 

( V XI* V YI» V zi) 

VXI, VYI, VZI 

inertial velocity vector 
and its magnitude, mps 
(fps) 

V I 


VELI 

magnitude of V , mps (fps) 

— IG 


U, V, W 

inertial velocity in the 
G-frame, mps (fps) 

-It 


VXIT, VYIT, VZIT 

velocity of target 
vehicle in ECI system, 
mps (fps) 

h 


VELA 

relative velocity, mps 
(fps) 

ho 


UR, VR, WR 

relative velocity in the 
G-frame, mps (fps) 

hi' 

( V WCI* V RYI» V RZl) 

VRXI , VRYI, VRZI 

relative velocity vector 
In the inertlel frame, 


aps (fps) 

II-v 





Internal 


Definition 


Math symbo l 


^VE 


*WI 



V WYI* 



V 

^WG 


V 

P 

V 

QO 

u 


'jett 


PC 


w"** 

\ 


w, 


PR 


w 


stg 

["»]• ["•] 

II -10 


Fortran rvabol 

VXVE, VTVE, VZVE velocity of vehicle in 
vernal equinox system, 
■pa (fps) 


VWXI, VWYI, VWZI wind velocity vector In 



the Inertial frame, mps 
(fpa) 

vw 

wind velocity, mps (fps) 

UW, W, ww 

wind velocity vector in 
the G-frame, ops (fpa) 

PGVEL 

perigee velocity, spa 
(fps) 

KYPVEL 

outgoing asymptote 
velocity, mps (fpe) 

WD0T 

total time rate of change 
of vehicle weight, N/s 
(lb/s) 

WEIC0N 

total weight of propel- 
lant consumed, N (lb) 

WEIGHT 

gross vehicle weight, N 
(lb) 

WJETTM 

jettison weight, N (lb) 

— 

weight of propellant con- 
sumed per phase, N (lb) 

— 

initial propellant weight 
N (lb) 

— 

maximum flowrata for the 
i*** engine, N/s (lb/s) 

WPR0P 

weight of propellant re- 
maining, N (lb) 

WCTSC 

vehicle stage weight, N 
(lb) 

Wl), W0PT, WE 

weighting matrices for u, 


!> and e 
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Math symbol 


Internal Fortran symbol Definition 


V y B’ *B 
“bR’ y BR* *BR 
x cg* y cg* *cg 


V y G* *G 


X I* y I» *1 


x 


1 


— coordinate axes of the 
body frame 

coordinate axes of the 

body reference frame 

XCG, YCG, ZCG coordinates of the center 

of gravity In the body 
reference system, a (ft) 

conponents of a vector In 

the geographic frame, m 
(ft) 

XI, YI , ZI conponents of the radius 

vector in the inertial 
frame, m (ft) 

— general state variable 


V y L* 


x 

~n 


coordinate axes of the 

launch frame 

GINTJ state vector at the n^ 

event 


V y R* 


z 


R 


x i*f’ y ref ’ *ref 


X. 

a, 8, o 


®T 


conponents of the radius 

vector in the planet 
frame, a (ft) 

XREF, YREP, ZREF coordinates of the aero- 

dynamic reference point 
in the body reference 
system, a (ft) 

DGENV general dependent variable 

ALPHA, BETA, aerodynamic angle of at- 

BNKAUG tack, sideslip, and bank, 

red (deg) 

ALPT#T total angle of attack, 

rad (deg) 
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Math symbol 

Internal Fortran symbol 

Definition 

\ 

\ 

Y I* y r* y a 

GAM4AI, GAMMAR, 
GAMMAA 

Inertial, relative, and 
atmospheric relative 
flight path angles, rad 
(deg) 

/ 

t 1 


GAMMA(I) 

step-slxe parameter on the 
j th trial step 


AE 

— 

increment in eccentric 
anomaly, rad (deg) 

> 

Ah 

— 

increment in altitude, m 
(ft) 

« 

At 

DT 

increment in time or Inte- 
gration step slse, s 


AV 

DV 

increment in velocity, 
mps (fps) 


* 

AV 

VIDEAL 

ideal velocity, mps (fps) 

t 

\ 

av a 

DLR 

atmospheric velocity loss, 
mps (fps) 

\ 

AV 

c 

DVCIR 

velocity required to 
circularise an orbit, 
mps (fps) 


av e 

av g 

DVEXS 

GLR 

excess velocity, ape (fps) 
gravity loss, aps (fps) 


av m 

DVMAR 

velocity margin, aps (fps) 


av p 

ATLR 

atmospheric pressure 
loss, aps (fps) 


av tv 

TVLR 

thrust vector velocity 
loss, mps (fps) 


4 a 

RTASC 

right ascension, rad (deg) 


6 « , u 

con* , dock 

SCONE, S CLOCK. 

cone and dock angles of 
Sun vector in body 
system, red (deg) 
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Internal Fortran symbol Definition 


n 

ETA 

engine throttling param- 
eter 

e 

L0NG 

planet relative longitude 
rad (deg) 

* 

e 

— 

longitude reference, rad 
(deg) 

e i 

L0NGI 

inertial longitude, rad 
(deg) 

V V ^ZL 

L0NL, LATL, A2L 

longitude, latitude, and 
azimuth of L-frame, rad 
(deg) 

0 

mix 

TRUNMX 

maximum true anomaly for 
hyperbolic orbit, rad 
(deg) 

8t i 

TRKLNI 

longitude of tracker 1, 
rad (deg) 

X 

A2REF 

azimuth reference, rad 
(deg) 

X 

STPMAX 

maximum admissible step 
size for the iteration 
algorithm 

u 

MU 

gravitational constant, 
a 3 /s 2 (ft 3 /e 2 ) 

V 

__ 

index 


i 
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Math symbol 


Internal Fortran symbol 


Definition 


Aggy argument of vehicle (l.e.» 

angular location of ve- 
hicle, measured from 
ascending node in orbital 
plane), rad (deg) 


p(h) 

— 

atmospheric density, kg/m 3 
(slug/ft 3 ) 

T 

— 

trajectory propagation 

♦c 

GCLAI 

geocentric latitude, rat 
(deg) 

♦« 

GDLAX 

geodetic latitude, rad 
(deg) 

♦l» 

R0I 1 , YAWI, PITI 

inertial roll, yaw, and 
pitch measured as positive 
rotations from the L-frame, 
rad (deg) 

V ®R* *R 

YAWR, PITR, R0LR 

relative yaw, pitch, and 
roll, measured in a posi- 
tive sense from the geo- 
graphic freme, rad (deg) 

n 

LAN 

longitude of ascending 
node, rad (deg) 

«p 

0MEGA 

angular rotation rate of 
planet about the polar 
axis, rad/s (dsg/s) 

«s 

RASGM 

right as cent ion of 
Greenwich meridian, 
rad (deg) 

fa) 

— 

argument of perigee, red 
(deg) 

« - (v <*y «,) 

R0LBB, PITBD, 
YAWBD 

Inertial angular velocity 
components about the body 
axis, rad/s (deg/s) 

& 

R0LBDD, PITBDD, 
YAWBDD 

inertial v.igular acceleration 
components about the body 
axis, rad/s 2 (deg/s 2 ) 


11-14 


Hath symbol 


Inf rnal Fortran symbol Definition 


< >A 

refers to atmosphere rela- 
tive variables 

(> cg 

refers to center of 
gravity 

< >1 

refers to inertial 
variables 

< >n 

refers to n** 1 event 

() P 

refers to thrust applica- 
tion 

< >R 

refers to Earth-relative 
variables 

* ^Ref 

refers to aerodynamic ref- 
erence point 


refers to sea-level condi- 
tions 

< >v.c 

refers to vacuum condi- 
tions 

< >W 

refers to wind relative 
variables 

( >* 

refers to state from which 
downrange and crossrange 
are referenced; refers to 
optimal conditions 

(_) 

denotes vector quantity 

( )' 

denotes transpose of a 
vector 

e 

( ) 

denotes total derivative 
with respect to time 
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COORDINATE SYSTEMS 


POST uses numerous coordinate systems to provide the neces- 
sary reference systems for calculating required and optional data. 
These coordinate systems and the key transformations are described 
below. 


Coordinate System Definitions 

Earth-centered inertial (ECI) axes |xj, yj, Zj) •- This sys- 
tem is an Earth-centered Cartesian system with Zj coincident 
with the North Pole, Xj coincident with the Greenwich Meridian 
at time zero and in the equatorial plane, and y^ completing a 

right-hand system. The translational equations of motion are 
solved in this system (.fig* HT-1)* 

Earth-centered rotating (ECR) axes |x^, y^» z r)' _ s y s_ 

tem is similar to the ECI system except that it rotates with the 
Earth so that x is always coincident with the Greenwich Merid- 
ian (fig. III-l) . 

Earth position coordinates ® ^ j • These are the fa- 

miliar latitude, longitude, and altitude designators. Latitude 
is positive in the Northern Hemisphere. Longitude is measured 
positive East of Greenwich. Altitude is measured positive above 
the surface of the planet (fig* III-l). 


Geographic (G) axes ^Xg, y G> * - This system is located 

at the surface of the planet at the vehicle's current geocentric 
latitude and longitude. The x G axis is in the local horizontal 

plane and points North, the y G axis is in the local horizontal 

plane and points East, and z G completes a right-hand system. 

This system is used to calculate parameters associated with azi- 
muth and elevation angles (fig* III— 2), 


Inertial launch (L) axes (x L , y L> * L )- - 


This is an iner- 


tial Cartesian system that is us7d as an inertial reference 
system from which the inertial attitude angles of the vehicle are 
measured. This coordinate system is automatically located at the 
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geodetic latitude and inertial longitude of the vehicle at the 
beginning of the simulation unless overridden by user input of 
LATL and L0NL. The azimuth, A^* is zero unless overidden by 

user input. The orientation of this system is such that is 

along the positive radius vector if $ is input as the geocen- 
trie latitude, or along the local vertical if ^ is not input 
or is input as the geodetic latitude, is in the local hori- 
zontal plane and is directed along the azimuth specified by A , 
and completes a right-hand system. This system is intended 

foi use in simulating ascent problems for launch vehicles that 
use either inertial platform or strapdovm-type angular commands. 
The inertial angles, 0^ are always measured with 

respect to this system and are automatically computed regard 1 ess 
of the steering option (1GUID) being used (fig, II1-2). 

Body (B) axes (xg, y R , z R j.- The body axes form a right- 

hand Cartesian system aligned with the axes of the vehicle and 
centered at the vehicle's center of gravity. The Xg axis is 

directed forward along the longitudinal axis of the vehicle, y 

B 

points right (out the right wing), and z D points downward, com- 

D 

pleting a right-hand system. All aerodynamic and thrust forces 
are calculated in the body system. These forces are then trans- 
formed to the inertial (I) system where they are combined with 
the gravitational forces (fig. III-3) 







j 
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Body reference (BR) axes (x ]JR , y fiR , z RR ) The body refer- 
ence system is a right-hand Cartesian system aligned with the 
body axes as follows. The x feR axis is directed along the nega- 
tive Xg axis, the y fiR axis is directed along tha positive 
y axis, and the z RP is directed along the negative z axis. 

This system is used to locate the vehicle* s center of gravity , 
aerodynamic reference point , and engine gimbal locations for the 
static trim operation (fig. HI-3) 


Orbital elements (h a> h p , i, u>, 0, a).- This is a nonrectangu- 

lar coordinate system used in describing orbital motion. The or- 
bital elements are apogee altitude, perigee altitude, inclination, 
longitude of the ascending node, true anomaly, and argument of 
perigee. The apogee and perigee altitudes replace the standard 
orbital elements of semimajor axis and eccentricity (fig. III-4) . 


Vernal Equinox <VE) Axes (x^, y^, ZyJ.- T** 18 is the 1950 

mean equator and equinox Earth centered inertial system. The Xy^ 

axis is in the equatorial plane and is directed forward of the 
vernal equinox of 1950, the axis is directed along the north 

pole, and y completes the right hand system (fig. I1I-1) . 

VJj 
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Attitude Angles 

The program contains the following standard attitude refe 
ence systems: 

1) Inertial Euler angles; 

2) Relative Euler angles; 

3) Aerodynamic angles; 

4) Inertial aerodynamic angles; 

These variables are defined and illustrated below: 

1) Inertial Euler angles (fig. III-5): 

♦i " Inertia l ro11 angle. The roll X B 

angle with respect to the L- I 

frame (first rotation), y L l\ ^ 

'Pj ~ Inertial yaw angle. The yaw ** ♦j* I y B ^ 

angle with respect to the L- y ' 

frame (second rotation) , ^ 

9j - Inertial pitch angle. The pitch 

angle with respect to the L- Figure III-5.- Inertial Euler Angles 
frame (third rotation) ; 

2) Relative Euler angles (fig. III-6): 

“ Relative yaw angle. This is *B 

the azimuth angle of the \ 

axis measured clockwise from x 
the reference direction (first G 
rotation) , 

0 R ~ Relative pitch angle. This is 
the elevation angle of the x^ 

axis above the local horizontal / '***-£*. y g 

plane (second rotation) , * 

z 

B 

* R “ Relative roll angle. This is 





the roll angle about the x_ v< ...... m * ... 

e Figure III-6.- Relative Euler Angles 

axis (third rotation). 



3) Aerodynamic angles (flgc III-7): 

o - Bank angle. Positive a is a 
positive rotation about the 
atmosphere relative velocity 
vector (first rotation) , 


0 - Sideslip. Positive S is a nose- 

left (negative) rotation when 
flying the vehicle upright (sec- . 

ond rotation) , y_ 

B 

a - Angle of attack. Positive a 

is a nose-up (positive) rotation 
when flying the vehicle upright 
(third rotation); Figure 111-7*- Aerodynamic Angles 

A) Inertial aerodynamic angles (fig* I I 1-8): 



o^ - Bank angle* Positive o ^ is a 

positive rotation about the 
atmosphere inertial velocity 
vector (first rotation). 


3 - Sideslip* 


Positive is a nose 


left (negative) rotation wh*n 
flying the vehicle upright (sec- 
ond rotation), 



*B 

I 


a 


I 


Angle of attack. Positive a 

1 Figure III-8*- Inertial Aerodynamic 
is a nose-up (positive) rotation Angles 

when flying the vehicle upright 
(third rotation); 
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Transformations 


Numerous matrix transformations are required to transform 
data between the coordinate systems described in the previous 
section. The most important of these transformations is the [IB] 
matrix. The inverse (transpose) of this matrix is used to trans- 
form accelerations in the body frame to the planet-centered in- 
ertial frame. The remaining transformations are generally used 
to either compute [IB] or to transform auxiliary data into some 
convenient output coordinate system. 

The [IB] matrix is functionally dependent on the attitude 
of the vehicle. This dependence is described by equations re- 
lated to the attitude steering option selected by the user. The 
following matrix equations, which depend on this steering option, 
are used to compute the [IB] matrix. 

[IB] * [LB] [IL ] (body rates or inertial Euler angles) \ 

[IB] = [GB] [IG ] (relative Euler angles) \ (III-l) 

[IB] * [AB][GA][IGj (aerodynamic angles) ) 

The basic relationships between the coordinate systems de- 
fined by these equations are illustrated in figure III-9. The in- 
verse transformation can generally be computed by merely trans- 
posing the matrix elements because of the orthonormality of 
these matrices. 



Figure UI-9'~ Matrix Transformations 
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A summary of these matrices is given below. The symbols s 
and c denote sin and cos, respectively. 


[IL] 


[IL], inertial to launch.- 

- The [IL] 

matrix depends 

on * L 

and A^, and is given by 




c *l c0 l 

c *l s9 l 


s4 l 

8 *L c9 L SA ZL - ^ZL^L 

^ZL^I. + 

sA ZL S V e i 

" sA ZL Cf 

' sA ZL S0 L ~ ^ZL^L^L 

sA zl c6 l ' 

“zl'V'i 

cA zl c * 
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[LB], launch to body .- The [LB] matrix is computed indi- 
rectly from the body rates by integrating the quaternion equa- 
tions, or directly from inertial Euler angles. When the body 
rate option is used, the quaternion rate equation 


eo 

*1 

*2 

*3 


1 

2 


~ e i 

e 2 

e 3 

eo 

e 2 

-e 3 

eo 

~ e l 

®3 

eo 

e l 

-e 2 


CHI-3) 


is integrated to compute the [LB] matrix, which is then given 
by 


[LB] - 


e 2 + e 2 - e 2 - - 2 
*0 + e l *2 1 

2(eje2 - eoe 3 ) 
2(e 3 e 3 + eoe 2 ) 


2(e 3 e 2 + eo e 3> 

- e 2 + e 2 - e 2 
2(02«3 " «0*l) 


2(e 3 e 3 - e 0 e 2 ) 
2(eoej + e 2 e 3 ) 

.2 . ,2 . ,2 i e 2 

e 0 *1 e 2 * e 3 
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When the inertial Euler angle option is used, [LB] is computed 
directly as 


[LB] = -s. 


C 'L C( I 


+ SiijSOj 

s v j si ^^c0 ^ 

- C^j-SCj 

- S -l 

c; l Cv l 


SVjCy 


C 'l S ’l 

C 1 S -1 S °I 

- S^CUj 

s$ I si fI se I 

+ 


(IH*5) 


[IG], inertial to geographic .- The [IG] matrix depends on 
the geocentric latitude and the inertial longitude, and is giver 
by 


[IG] * I -s« 


-S<J> C0 T 
c I 


— s s0 T 
c I 


(III-6) 


-c4> c T 0 -c* s6 t -s$ 
cl cl c 


[GB], geographic to body .- The [GB] matrix depends on the 
relative Euler angles, and is given by 


C V*R 


[GB] « ^R^R^R “ ^R^R 


C V*R 


S$ R S0 R SlJ> R + C$ R C^ R 


C V e R c *R + 8 V*R c W^R - 8 V*R 


S V e R (IIW) 


C V e R 


[GA], geographic to atmospheric relative velocity system 
(ARVS) The [GA] matrix depends on the atmospheric relative 
flight azimuth and flightpath angles, and is given by 


cY a cX a 


[GA] - -s\. 


*V X A 


C W 


*v\ 
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JAB] > ARVS to body ,- The [AB] matrix depends on the aerody- 
namic angles, and is given by 


cac3 

[AB] » s8 


-casSc^ + sotsa 
c8co 

-sasBca - cas a 


-cotsBsa - sotco 
cBsa 

-sasSso + catca 


(in-9) 


Other transformations, which are not related to the calcula- 
tion of the [IB] matrix, are presented below. 

[IP], inertial to planet relative .- The [IP] matrix trans- 
forms between the Earth-centered inertial frame and the Earth- 
centered rotating frame. This matrix depends on the rotation 
rate of the planet and the total elapsed time of flight, and is 
given by 


eft t 

s ft t 

0 


P 

P 


H 

-sft t 
P 

eft t 
p 

0 

(in-io) 

. 0 

0 

1_ 

• ! 


[RB], body reference to body .- The [RB] matrix transforms 
data in the body reference system to the body frame. This matrix 
has a constant value and is given by 


(III-ll) 


between 

the ECI 

frame 

and 

the vernal equinox 

by 

m 






c 6 

8 

0 

0 


[IV] - 

-s 0 

c 

0 

0 



_0 

0 


1^ 

» 

where 

0 - CHA 

♦ ft 

(*- 

t ref erence ) * 


( 111 - 12 ) 


III-10 


IV. PLANET MODEL 


The planet model Is composed of three types of data and equa- 
tions. These are: (1) oblate planet geometry and constants, (2) 

an atmosphere model that computes atmospheric pressure, density, 
temperature, and speed of sound, and (3) a gravitational model 
that computes the gravitational accelerations. The user selects 
the appropriate models and Inputs the corresponding data. The 
input data and the equations used in these models are described 
below. 


Oblate Spheroid 

The 1960 Fisher Earth model is preloaded into the program. 
This model is defined by the equatorial radius Rg, the polar 

radius Rp, the rotation rate the gravitational constant 

u, and the second, third, and fourth gravitational harmonics, 
J 2 » J 3 . and J 4 , respectively. The stored values for these 

constants are: 


Rg - 2.0925741 x 10 7 ft, 

Rp - 2.0855590 x 10 7 ft, 
fl p - 7.29211 x 10“ 5 rad/s, 

P - 1.40/ 6539 x 10 16 ft 3 /s 2 , 
J 2 - 1.0823 x 10- 3 , 

J 3 - 0, 

J 4 - 0. 


The constants J 3 and J 4 are preloaded as zero, but can be ini- 
tialized by input. For example, if the Smithsonian Earth model 
is desired, then these constants would be input as 


J 2 - 1.082639 x 10“ 3 , 

J 3 - -2.565 x 10“ 6 , 

J 4 “ -1.608 x 10" 6 , 
u - 1.407645794 x 10 16 ft 3 /s 2 . 
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n p - 7.29211515 x 1<T 5 rad/s, 

Rg - 2.092566273 x 10 7 ft, 

Rp - 2.08550242 x 10 7 ft. 

The geometry of this spheroid is illustrated in figure IV-1. 
The pertinent equations related to this model are 

♦c • • 1 “"‘ (*iAi) 

■ tan' 1 |k tan $ c ), k « 

R. - «E (l + (h - l>.in 2 ♦ c )-’ s 

h ‘ r l ' V 

where 4> is the geocentric latitude, 4> is the geodetic lati- 
c 8 
tude, 6 p is the inertial longitude, 0 is the relative longi- 
tude with respect to the planet, Tj is the distance from the 
center of the planet to the vehicle, R g is the distance from 

the center of the planet to the planet surface, and h is the 
distance from the planet surface to the vehicle. 



North pole 



• Vehicle 


South pole 

Figure IV-1.- Oblate Planet 


Gravitational Model 

The gravitational model includes optionally second, third, 
and fourth harmonic terms. The potential function for this model 

is 



1 


(IV-?) 
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Th« gravitational accelerations calculated from this potential 
function are: 

c .-12- 

TU 3x, 


-w P (*, r) 


^YI 3y, 


-y ^ P (a, r) 


G - - 12 - 

ZI 3z, 


- - p- (l + JR 2 (3 - 5Z 2 )| z + H ~~ ^6z 2 - 7z 2 Z 2 - | r 2 J 

+ DR 4 | — - iOZ 2 + 9Z 4 j zj , 


where x - Xj, y - y Jt z - a , r - r , and 


R " hn 


Z - * T /r. 


J - f J 2 


H - | J 3 


D - - f Ji, 


1 + JR 2 (1 - 5Z 2 ) ♦ H ~ (3 - 7 Z 2 )a 


+ DR 4 (9Z 4 - 6Z 2 + 


L 
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Atmosphere Models 


POST has the optional capability of three atmospheric models — 
the general table lookup, the 1962 U.S. standard atmosphere, and 
the 1963 Patrick AFB atmosphere using polynomials. The general 
table lookup model gives the user the flexibility of inputing his 
own atmospheric model if none of the preloaded models is adequate. 
This is particularly useful in performing trajectory analysis for 
planets other than Earth. The parameters required to define the 
atmospheric effects are the atmospheric pressure p, atmospheric 
density o, speed of sound C^, and atmospheric temperature 

T. These parameters are functions of the oblate altitude h. 

Table lookup atmosphere model .- The table lookup atmosphere 
model can be defined entirely by using tables that show pressure, 
temperature, speed of sound, and density as functions of altitude. 
The speed of sound and density tables can be omitted if desired; 
in this case, the speed of sound and density are computed as 


c s 


K £ 
*•2 j » 


(IV-5) 


where 


K, 

Kl M- 


R* 


y » ratio of specific heats 
Mq * molecular weight 

R* * universal gas constant . 

196.2 U.S. standard atmosphere modal .- The 1962 U.S. stand- 
ard atmosphere model la given as a function of geopotential alti- 
tude J, which is computed as 


R A h 
R A + h* 


(IV-6) 
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F 


7 


where 


k - average Earth radius * ^ + ^p) 

A 

h - oblate altitude. 

The molecular scale temperature, T^, is defined by a series 
of linear segments as a function of geopotential altitude 


(“»)• 


The corner points connecting the straight-line segments aie 
referred to as base altitudes (Hg), base temperatures jT^j, 

etc. From a table of base altitudes, base temperatures, and 
dT )dH /L^ (the slope within the linear segments), the tempera- 
ture at any desired altitude can be calculated from the following 
equation: 


'M - \ + \ ("g ' V 


Values of P B> T , *md 1,^ versus are presented In 


table IV-1. 

The atmospheric pressure is determined as follows: 



’ T 


( 8 o M o \L 

1 

p - P B 

T h 

exp 

\ p* f 



L j 

r 

i 

- 

J 


P - P B exp 


«o M 0 1 “- h b) 


R* 




for segments with Lj^ 0, ® n< ^ 

for segments with “0, 


where P^ is the bese pressure corresponding to the given base 

altitude H . These base pressures can be calculated once the 

8 


sea- level pressure, P Q , and the temperature profile have been 
specified. 


Having calculated the temperature and pressure, 
p, speed of sound, C g , and atmospheric viscosity, 

determined as follows: 



the density, 

V are 


(IV-9) 


where is the acceleration of gravity at sea level, is 

the molecular weight of air at sea level, R* is the gas con- 
stant, y is the ratio of specific heats, and 8 and S are 
Sutherland's constants. 


M q - 28.9644 


R* 


Y 

B 


8.31432 • 10 
1.40 

1.458 x l(T fc 


3 (°K) 


_J 

(kg-mol) 


k g 

sec m (°K) % 


) 


\ 
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S - 110. 4°K -* 198 . 72® R 
g Q - 9.80665 m/sec- - 32.174 ft/sec . 


In the 1962 U.S. standard atmosphere, the molecular weight 
varies with altitude above approximately 90 km; in POST the molec- 
ular weight is assumed constant, resulting in a slight discrepancy 
above 90 km. in the 1962 U.S. standard atmosphere, geometric alti- 
tude is transformed to geopotential altitude, which used through- 
out. Thus, above 90 km, a constant slope of molecular scale tem- 
perature versus geopotential altitude is used instead of the con- 
stant slope of temperature versus geometric altitude. 
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Table IV- 1. - 1962 U. S. Standard Atmosphere Profile 


h b* ft 

P B» P sf 

V OR 

* R/ft 
B 

0.0 

C. 21162166 + 4 

518.67 

-0.35661600 - 2 j 

36 089.239 

0.47268050 + 3 

389.97 

0.0 

65 616.797 

0.11434543 + 3 

389.97 

0.54863995 - 3 

104 986.87 

0.18128948 + 2 

411.57 

0.15361920 - 2 

154 199.48 

0.23163263 + 1 

487.17 

0.0 

170 603.68 

0.12322603 + 1 

487.17 

-0.10972801 - 2 

200 131.23 

0.38032532 + 0 

454.77 

-0.21945600 - 2 

259 186.35 

0.21673064 - 3 

325.17 

0.0 

291 151.57 

0.34333824 - 2 

325.17 

0.16953850 - 2 

323 002.74 

0.62814785 - 3 

379.17 

0.28345707 - 2 

354 753.59 

0.15361733 - 3 

469.17 

0.56867005 - 2 

396 406.39 

0.52676024 - 4 

649.17 

0.11443751 - 1 

480 781.04 

0.10566108 - 4 

1 729.17 

0.86358208 - 2 

512 046.16 

0.77263469 - 5 

1 999.17 

0.57749093 - 2 

543 215.48 

0.58405376 - 5 

2 179.17 

0.40610461 - 2 

605 268.45 

0.35246030 - 5 

2 431.17 

0.29274135 - 2 

728 243.91 

0.14559124 - 5 

2 791.17 

0.23812804 - 2 

939 894.74 

0.39418091 - 6 

3 295.17 

0.20152600 - 2 

1 234 645.7 

0.84380249 - 7 

3 889.17 

0.16354849 - 2 

1 520 799.4 

0.22945543 - 7 

4 357.17 

0.11010085 - 2 

1 798 726.4 

0.72259271 - 8 

4 663.17 

0.73319725 - 3 

2 068 776.5 

0.24958752 - 8 

4 861.17 

0.0 


IV-8 



1963 Patrick AFB atmosphere using polynomials .- In this 
model, pressure and temperature are calculated as functions of 
geometric altitude (h) . These parameters are calculated in met- 
ric units and converted to English units if required. 

Pressure: 

1) Altitude region * 0 to 28 000 meters : 

P - P! exp (A + A x h + A 2 h 2 + A 3 h 3 + A 4 h 4 + A 5 h 5 ) 
where P} * 10.0 Newtons/cm z ; 

2) Altitude region * 28 000 to 83 004 meters: 


3) Altitude region * 83 004 to 90 000 meters: 


P B axp 


-1.373301523 * 10 12 h - h \ 
T b (6344860 + h) (6344860 + h B )j ; 


4) Altitude region = 90 000 to 700 000 meters: 


L n (P ^ “ L n < P b) + L (6344860 

m 


1.373301523 * 10 12 \ 

44860 + h) v 6344860 + h ) 1 


T 

(\ + L » ( h - 


1V-9 




Temperature: 

X) Altitude region * 0 to 10 832.1 meters: 

T * T* * A + Ai h + A2 h 2 + A3 h 3 + Ai+ h 4 + A5 h 5 ; 

2) Altitude region * 10 832.1 to 83 004 meters:* 

T * A + Ai h + A 2 h 2 + A 3 h 3 + A 4 h 4 ■+ A 5 h s ; 

3) Altitude region * 83 004 to 90 000 meters: 

1 ■ T B + he ( h - Ns)' 

However, in this region L k ” and thus 

T ■ T * 180. 65°K; 

B 

4) Altitude region ■ 90 000 to 700 000 meters. 
T “ T M“ T M B + L m( h " h B)- 


> (1V-12) 

i 


Density: 

1) Altitude region ■ 0 to 28 000 meters : 

p - Pl exp (A + Aj h + A 2 h 2 + A 3 h 3 + A 4 h" + A 5 h b ) ; 

2) Altitude region - 28 000 to 700 000 meters: 


(1V-13) 


p ■ (34.83676) 


^Virtual temperature is the seme ss kinetic temperature 
above the 10 832.1-meter altitude. 
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Derived Coefficients for the 1963 Patrick AFB Atmosphere Model 
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Table IV-3. - 1963 Patrick AFb Molecular 
Temperature Profile and Gradient Profile 


V km* 


90 

100 

110 

120 

150 

160 

170 

190 

230 

300 

400 

500 

600 

700 


Tm OK 

V 


180.65 

210.65 

260.65 

360.65 

960.65 
1 110.65 
1 210.65 
1 350.65 
1 550.65 

1 830.65 

2 160.65 
2 420.65 
2 590.65 
2 700.65 


V CK 


3.C 

5.0 
10.0 
20.0 
15.0 
10.0 

7.0 

5.0 

4.0 
3.3 
2.b 
1 . 


l.: 


*Altitude range: 90 000 to 700 000 meters 



Pressure and density ratios: 


Altitude region * 0 to 700 000 meters: 


D 



Velocity of sound: 

V„ - (20.046707) (T)’ 2 . 

w> 


(IV-14) 


(IV-15) 


The atmosphere model-derived coefficients are presented in 
table LV-2. The molecular temperature gradient is documented in 
table LV-3 for geometric altitudes from 90 to 700 km. 


Winds 


The atmospheric wind velocity components are input in tables 
using either meteorological or vector notation. If these tables, 
which are normally functions of oblate altitude, ere not input, 
then the atmosphere is assumed to rotate uniformly with the 
planet . 


The wind velocity components can be input directly in the 
geographic frame by defining Uy, v w , and or by defining the 

wind speed ( v w ). the vertical component ( w w j » the wind a zimuth 
( A zw) * and the wind azimuth bias ( A zWb) * resulting wind velo- 

city components in the G- frame are: 


^wc 


V w <h) cos (A zw (h) + A zwb ) 

v w (h) 8in ( A zw (h) + a zwb) 


(1V-16) 


L“w (h) 


IV-14 


V * V - j'l * r - V. 

— AI -I ~P “I -WI 


and its magnitude is given by 


V A = j4l * — AI 


(IV-18) 


(IV-19) 


It is clear from the abo^e equation that in order to input vector 

wind data A„,, n must be input as zero, whereas for meteorologic 
ZWB 

data the preloaded value of 180° should be used. 

The wind velocity in the EC1 frame is then given by 

V WI * W' -WG 

Thus, the atmospheric relative velocity vector in the ECI frame 


(1V-17) 
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V. 


VEHICLE MODEL 


The various physical properties of the vehicle are modeled 
by the user when he selects the pertinent options from the set of 
vehicle simulation modules. The equations used in these modules 
are presented below. 


Mass Properties Model 


The gross weight of the vehicle at the beginning of each 
phase is given by 

Vf = w + w 

G stg pld 


(V-l) 


where M is gross weight without payload and 
stg 


pld 


is the pay- 


load weight. For phases other than the first, the gross weight 
can optionally be computed as 


“g ■ “c - H )eu - W PR- 


(V-2) 


where W is the gross weight on the positive side of the cur- 

G 

rent event, W is the gross weight on the negative side of the 

G 

current event, W. _ is the jettison weight, and W is the 
jett 

weight of propellant remaining. These options are obtained au- 
tomatically, based on user input. 

The propellant remaining is given by 

VL 


PR 


\ - W pc ’ 


( V— 3 ) 


where is the initial weight of propellant and W is the 

1 

amount of propellant consumed. This latter term is given by 


w pc ■ /“ 


dt + W, 


( V-A) 


where W is the total rate of change of the vehicle's weight 


* 







V-l 


T 


At the beginning of each phase, the constant W^, can be 

either input or carried across the event as the total amount of 
weight consumed in the previous phase. 

The amount of propellant jettisoned can be computed as: 

1) The amount of propellant remaining at the beginning of 
the current phase, 

2) The amount of propellant remaining at some prescribed 
prior event. 

Tne constant jettison weight is either computed from an in- 
put constant value or determined from an input mass-fraction table. 
When a mass— fraction table is used the jettison weight is given 
by 


w. „ * w. 7-1 , 
jett P. [A J 


wnere A is the mass fraction computed from the table. 


Propulsion Calculations 

POST can simulate both rocket and Jet engines. Tne program 
can simulate up to 15 engines In either mode. 

Rocket engines.- There are two input options for engine data 
in the rocket mode. In the first option, tables for vacuum thrust 
and maximum weignt flowrate are input for each engine. In the 
second option, tables for vacuum thrust are input, along with 
tne vacuum specific impulse for each engine. The vacuum speci 
fic impulse is then used to calculate the mass flowrate. 

The rocket thrust per engine is given by 
\ ■ " T v.c t - \ p, " ) ' 


‘At Am, 


where n i s the throttle setting, T vac is the. vacuum thrust 

of the i th engine, Aj. is the nozzle exit area, and p(h) is the 

atmospheric pressure. Summing over all engines yields the total 
rocket thrust 

H 


eng 

L V 

i=l 


(V-7) 


where N is the number of thrusting engines, and N - 15. 

eng B 

The weight flowrate in the rocket mode is given by 

N 


W = 


eng 


E («r 

.) 

i-i 


N 


eng 


( TV3C / 

/ Is Pvac)i 

i-1 \ / 


jet engine 

mode the net jet thrust per 


(V-3) 


engine is given by 


(r) ■ ;(M - ”>• 


(V-9) 


where 


p(h) /l 


SL 


(V-10) 


V-3 


a monovariant table. The total jet thrust is then 


and (M, 

Cj 

given by 


n «« 

Tj - n 5^P( h >/(P SL ) (Tj/M- 

i«l 

The weight flowrate in the jet engine mode is 
N 

eng 


(V-ll) 


(V-12) 


The thrust vector components for both rocket and jet engines 
are determined from the thrust magnitude T or T, and the 

R i J i 

thrust incidence angles X and i . The thrust accelerations 

p i y i 

in the body axes are then given by 


eng 


-T3 




i*l 


cos i cos 1 
y-t 


sin i 


cos i sin i 
y i p i 


h 

m 


(V— 13) 


In the abcve equation T^ denotes the thrust magnitude for the 

i engine (either rocket or jet) and m denotes the instantaneous 
mass of the vehicle. The engine gimbal angles are determined 
from the static trim equations in the moment balance option or 
by input if the moment balance option is not used. The engine 
gimbal angles are illustrated in figure V-l. 






1 


r 

\ 





l 


r 


5i 



Figure V-l.- Engine Gimbal Angles 

Note that thrust misalignments can be simulated by inputting the 
engine gimbal angles and using the standard three-degree-of- 
freedom option. 


Aerodynamic Calculations 


The aerodynamic force coefficients can be expressed in terms 
of the lif*-, drag, and side-force coefficients C L » C^, and C . 

(fig. V-2), where C L and Cj, are directed normal to, and along the 
velocity projection in. the x^,-Zg plane. Note that produces 

a side-force, F , acting in the direction of y fi . 


L.ift and drag force coefficients are transformed to axial 
and normal force coefficients as follows: 


V 

3 , 

cos a 

-sin a 

C N 


s Ln a 

cos a 




- 


whera a is the angle-of-attack. 



(V-14) 


V-j 



Figure V-2.- Aerodynamic Angles 

The aerodynamic coefficients can also be expressed in terms 
of the axial force, normal force, and side force, C , C. , and C V9 

A iS I 

respectively. Here and produce forces that act in the 
-* B anc * directions, and produces a force acting along y^. 

Each aerodynamic coefficient is computed by interpolating 
the values in the table. In general, eight tables are allocated 
to each coefficient. These tables can be monovariant, bivariant, 
or trivariant, and seven tables per coefficient can have arbi- 
trary Hollerith mnemonic multipliers. This generality enables 
all standard forms of aerodynamic data to be directly input into 
tue program. 


V-6 


The aerodynamic force coefficients are given by: 


S = V + C D M C, + S. \ +C D. M C_ 

0 D : P D 5p dy d 5y 


C, “ C L + C L M c + C M 
“ L 0 L 0 C L 0 L 5p C L ip 


C A * V + C A \ + V \ * \ \ 

0 A 6 p A 5 5y 


(V- 13) 


C„ - C._ + C M + C M 

W C XT W- c 

0 N o A , 

P - 


c * ■ S ♦ S % * S s % 

y °, 


the aerodynamic moment coefficients are given by: 


C M - C. f + C.. M„ + C M 

rl il n *i ow M. C 

0 M op M,. 


C - C + C M + C M 
n n C M n. C 

0 .i Oy n . 


(V-16) 


In the above expressions, C , Cp, v C Q , etc, denote 

U( ) 6p 5y 

the tables, and M-, , M , M , M , etc, denote the mnemonic 

h)„ °D D r U x 

0 5p 5y 

table multipliers. Typical table arguments and multipliers would 
be a, 0, li, R.^, 6p, and 6y. 

The liach number and dynamic pressure are given by: 

H - V A 1 


1 ‘ 2 0 V A’ 


(V— 17) 


where p is the atmospheric density, V A is the velocity of the 
vehicle with respect to the atmosphere, and C g is the speed of 

sound, These atmospheric parameters are determined from the at- 
mospheric models as a function of the altitude h above the oblate 
spheroid; 1. e. , 

V- 7 


p n y * ' ,' wp pii ii l u ii m m mi 



(V-22) 


The resulting accelerations in the body system are then ob- 
tained from 


Km 


— F 
m ^AB 


Aeroheating Calculations 

POST provides for a wide variety of aeroheating calculations. 
Some of these options are specific in nature and apply only to 
particular vehicles, whereas others are quite general. The gen- 
eral heat rate option is based on trivariant table interpolation 
and provides complete flexibility with regards to vehicle shape 
and heat-transfer methodology. The various heat rate equations 
are described below. 


W 

■% 

■5 


'i 


i 

| 


} 


i 

i 
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Heat rate equations.- 


1 ) 


Chapman *8 equations . In this calculation the heat rate 
is given by 



where is the nose radius, p is the atmospheric density, 

and V is the reference circular orbital velocity, 
c 

2) General table lookup . This heat rate is given by 
Q - Q t (*l» *2. *3) . 


(V-23) 


(V-24) 


where x lt x 2 , and x 3 can be any internally computed variables. 
For example, the value® that would normally be selected are x\ - a, 
x 2 - h, and x 3 - V R> 

3 ) Modified Chapman's equation . Here the heat rate is given 
by 

Q - Q t (*i» * 2 . *3) Q c * 

where Q t is an arbitary table and Q c is the standard Chapman’s 
equation. 

4) Turbulent-flow heat rate . The turbulent- flow heat rate 
is given by 


Q 


Q t <*l» x 2» *3> 



(V-26) 


5) 


Maxim «« c<mterline heating . The equations for this method 
are ‘given ^elow in sequence . 

a) Altitude-velocity correction: 


Ah - 10 5 ^1.06112 - 6.16586 V^O^J 

+ 51.12090 (yi0‘\' - 20.66250 (v A/ /lO-) 5 j 

+ 22.52598 [V/J 10 '*) 2 ~ *8.2808<>| V A ^0 4 j 3 




) (V-27) 


J 


b) Maximum centerline heat rate at reference conditions: 

— if h . > 103 600 m: 

ref — 

«ref " 10 2 [277.93332 + 134.55760 h re£ /10 5 - 807.75941 (h. tf /10 5 ) 2 
+ 2.90536 (h re£ /l 0 s ) 3 + 722.36896 (l» re£ /10 5 ) 4 - 311.40176 

„ , „ (W 10 ’) I ; 

— if h , < 103 600 m; ' 

ref — 

q - 10 4 [ 7115 . 39692 - 34 881.13588 \ e{ /^ 5 + 69 844.23141 

— 71 534.98453 (h M1 /10 5 ) 3 + 37 506.13054 (>> r>f /10 5 ) ' - 

( h r.«/ l ° 5 ) 5 ]- 


- 8048.55112 


c) Angle of attack correction: 

q /q eno ■ [In (*)] 2 » 

n max , af max , a- 50 

where 

x - 10 2 [0.01136 + 0.01343 a/10 2 + 1.42672 (a/10 2 ) 4 - 0.75623 

(a/10 2 ) 5 ] + 0.30535 (a/10 2 ) 2 - 1.06269 (a/10 2 ) 3 . 

d) Maximum centerline heat rate: 


> CV-28) 


(V-29) 


Snax " (Saax.al/^lSnaxia-SO*) (^ref)* 

In addition to the heat rate calculations, the program also 
provides the capability to calculate other aeroheatlng indicators 
that can be used for trajectory shaping purposes. 

Aerodynamic heating indicators .- The heating rate for zero 
total angle of attack a^ is 


(V-30) 


Q - q V 


A* 


(V-31) 


V-ll 




The aerodynamic heating Indicator for zero total angle of attack 
is 


t 



0 


The heating indicator for non-zero angles of attack is 
given by 


Q' 


t 



9 


where 


and 


f M) 


/ 7 \ 5/7 

- I 1 + jM 2 sin 2 a' j K, 

| i + ^[ i -( i + 1 m2 * 1 " 2 a ') 2/7 ] 


K - 




% 


% 


% 




for a < 0° 


for a > 0* 


- Q' 

- 0 


- Q 

- 0 


■) 


. Q' 


for 8 < 0* 


for 8 > 0*. 
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The heating indicator for laminar flow Is calculated as 


K = f 17 600 K /— \ / 
iam J a T\ P 0/^ 


_k ) 

26 000 


where 


* * f (O. 

a T T 

The heating Indicator for turbulent flow is calculated as 


t0.8 / V 


1500 K 


°o) ( To“ooo j dt 


Ten-Panel Vehicle Heating Model .- Special aeroheating calct 
lations are available for a ten-panel vehicle model. The heatng 
ratios are referenced to the heat rate calculation. Tha total 
heat for each panel is given by 

9 i ■ H R 1 9 - 

where Q Is the total heat and H is the heat ratio for panel 

i 

i* The weight for each panel is the product of the weight per 
unit area and the area of the panel. The total weight is the sun 
of the individual weights for each panel: 


AW 

£ X 


where is the weight per unit area and as the area 

of the i 1 * 1 panel. 
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Steering Model 

The steering options control the attitude of the vehicle 
during the trajectory simulation. The general types of steering 
options available ares 

1) Body rates; 

2) Aerodynamic angles; 

3) Inertial Euler angles; 

4) Relative Euler angles. 

The body rates are generally used to simulate strapdovn-type 
systems and are computed from user-specified rate polynominals. 

K Srody»Mic ..gi.. at. ....rally u..a lor r..-try 

and the inertial and relative Euler angles are us tally used to 
simulate vehicles that employ inertial or local n "**“*“* ve *~ 
erence systems. All of these angles can be computed from. (1) 
polynominals; (2) tables; (3) piecewise linear functions; or (4) 
closed-loop linear feedback systems. 

The functional relationship used to compute the steering com- 
mands suggest two natural steering classifications: 

1) Rate steering; 

2) Angular steering. 

These classifications provide an efficient outline for pre- 
senting the equations used to compute the steering commands. 

Rate steering.- Rate steering uses the body rates, in con- 
junction with the quaternion equations, to determine the attitude 
of the vehicle. When using this option the user must specify: 

1) The initial attitude of the vehicle; 

2) The polynominal option used to compute the body rates. 

The Initial attitude is used internally to initialize the quater- 
nion rate equation 


i - i [E] u 
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Figure V-3*- Body Rates 

There are two options available for initializing the quaternion 
elements: (1) inertial Euler angles > and (2) aerodynamic angles. 

When inertial Euler angles are input the initial quaternion vector 
is given by 


•o - « (*!) * S. (*x) * • (»!> 
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whsre the asterisk denotes quaternion multiplication and where 


CO 

- COS 

(0.5 

♦«) 

+ sin 

(0.5 

♦0 
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(♦0 

- cos 

(0.5 

*0 

+ sin 

(0.5 

*0 

k 
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When aerodynamic angles are input , then the initial quater- 
nion vec*-*r is given by 

•o ■ 2 (* 2L ) * «W0) . . (♦,_) * . * . (9j) . . (-* c ) * .<-90) 

e (A^j * e (y^j * e(4) * e(-0) * •(<*)» < 


where 


s. ( a zl) “ 008 (°* 5 a zl) “ 8ln (°* 5 a zl) k 


_e(90) - cos (45) 


+ sin (45) 


e ■ cos (0.5 + sin (0.5 $ L ) j 

e (-0 L j - cos (0.5 6 L ) - sin (0.5 ej k 

ii (®j) " cos (®»5 + sin (o.3 k 

b ■ cos ( 0.5 - sin (0.5 $ c | j 

a 

%(a) - cos(0.5oj<4- (sin 0.5 o) i 

^(-0) ■ cos(o.50j- (sin 0.50) k 
^(o) " cos(0.5a|* (sin 0.5a) j. 


The user must also select the option for computing the body 
rates. These options are combinations of the basic rate poly- 
nominals shown at the top of next page: 


MX 





where a^, b^* and c^ are the polynomial coefficients* and 

x, y, and s are the polynomial arguments. 

The available combinations of these basic rate polynomials 

are: 

1) Input the coefficients and the arguments of 

2) Input the coefficients and the arguments of a, 6, a, 

and calculate id , and via 
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4> x *1 + d lO + (sin a) 6 

o)^ » &2 + d2^ + a 

I* 3 + d 3 a * ® ; 


3) Input the coefficients of the a, u , and u> poly 

noelals and calculate u> via ^ 7 8 

x 

w_ * ai + djd + tan a | 83 + djd - 

4) Input the coefficients of the a, and poly 

nomials and calculate u via 

y 

t w - at + (« - a 3 ) tan a 1 
dj + d 3 tan a J 

5) Input the coefficients of the u , u , and 6 poly- 

^ y 

nomials «»nd calculate u via 

z 

. d 3 r .1 

i *» 2 * a 3 - (cos a) 0 + — I - a* -(sin a) 6 I ; 

• m 

6) Input the coefficients of the a, a, and o> poly 

2 

nomials and calculate u> and u> via 

* y 

<»> x ■ a^ + d^d + tan a ^a 3 + d 3 d - ^ 

u ■ 82 + djd + a ; 


7) Input the coefficients of the a, and 0 poly- 

nomials and calculate u and u via 



. d 2 v -I 

“ y " a 2 + « + ^7 (« x - *1 " (sin o) B/ I 

d 3 , j 

a 3 + /«* x “ a l (sin a) ; J 


u * a» - cos 
z 3 
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8) Input the coefficients of the a, and 8 polynomials 


and calculate u and u via 
x. z 


■ aj + d]0 + (sin a) 8 
w ■ a 3 + d 3 o - (cos a) &; 


} 
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9) Input the coefficients of the 0 R , <J> R polynomials and 

and calculate w , u , and co via 
x y* z 


” — 


“ — 

U> 

X 


♦r " 8in 8 r *r 

Ui 

y 

* a + 

cos <J> R 0 r + sin cos 9 R i|> R 



cc* (j) R cos e R ij/ R - sin <j> R e R 
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The variables used in these equations are 






and 


u A « a + 
A 


v A - w + 
A 


du 

W Vi 

dh h 
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dv 

f) (r cos <p + r <t> sin c ) 
P\ c c c/ 

+ dh~ h 

dw 


— — h 


dh 



^[“a 4 a + Va + ”a“a] 


~w - 


dR 

— ♦ 
d<b T c 
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dR 
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d* ’ 
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sin <j> cos 4> 
c c 




Angular steering .- Angular steering uses four different 
functional relationships to determine the attitude of the vehicle. 
These functional steering equations are: 

1) Cubic polynomials ; 

2) Tables; 

3) Piecewise linear equations; 

4) Closed-loop linear feedback systems. 

Each of these steering techniques is available for all of the 
steering angles; i.e., aerodynamic angles, inertial Euler angles, 
and relative Euler angles. There is also a separate channel 
steering capability. This option enables the user to specify 
different functional relationships in each triplet of steering 
angles. This means, for example, that the an? le of attack could 
be computed from a polynomial, the bank angle from a table, and 
the sideslip angle by yet another method. However, this option 
does not allow the user lo mix the steering triplets (i.e., mix 
the aerodynamic angles with the Euler angles) . 


! 
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In the following discussion, it is convenient to let 0 de- 
note an arbitrary steering angle. The steering equations de- 
scribed in terms of this variable then apply to all of the steer 
ing angles. 

Polynomial steering: Under this option the steering angles 

are computed from a cubic polynomial 


£ V 1 ' 


where the user selects the coefficient and the independent 

variable y. The highest-order coefficient that is input deter- 
mines the degree of the polynomial. The argument can be selec- 
ted as any internally computed variable; e.g., time, velocity, 
altitude, etc. The constant term of the polynomial, c^, can 

be either input at the beginning of each phase or carried across 
as the value of the angle at the end of the previous phase. The 
polynomial coefficients are generally used as the independent 
variables for targeting/optimization. 

Table steering: Under this option the steering angles are 

computed via table interpolation, which is denoted by 

0(*) - 0 m T n [f(*)]. 


The user initially inputs the table multiplier 0 , 

m 


the order or 


interpolation n, and the table data (^, .f(y) )• The table mul- 
tiplier or the dependent table values can be used as independent 
variables for targeting or optimization. The order of interpo- 
lation can either be linear or cubic. The tables can be mono- 
variant, bivariant, or trlvarlant functions of the table argu- 
ments. 


Piecewise linear steering: Under this option the steering 
angles are computed from a general piecewise linear function of 
the form 


0(y) = c x 
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where c A is equal to 0 at the beginning of the current phase, 
c 2 is the desired value of 0 at a designated later event, y is 
equal to the value of the designated event criterion at the begin- 
ning of the current phase, y 2 is the desired value of the desig- 
nated event criterion at the designated event, and y is the cur- 
rent value of the designated event criterion. 

This option is similar to the polynomial option except that 
the values of 0 are specified directly rather than as 6 0 , 9, 

•• • • • 

9 and 0. Clearly, 0 is linear in time if y * t; otherwise 

0 is only linear in y. When the desired values of the steering 
angles are used as independent variables, the problem of cascaded 
steering effects is avoided and the targeting/optimization algorithm 
generally converges faster. This option also automatically com- 
putes the steering angle rates required to change the attitude 
to the desired value at the designated event, which reduces the 
problems related to guessing accurate initial pitch rates. 

Linear feedback steering: Under this option the steering angles 

are computed from the linear error-error rate feedback control 
law 


9 "=1 + S O'. - F d) + K R l r a - F d> < V - 65 > 

where cj is a nominal steering angle profile, is the dis- 

placement error gain, K is the rate gain, F - F , is the error 

R a q 

in the steering function F, and F' - D' is the derivative 
of the steering function error. 

This option i9, of course, the classic path control law, and 
enables the user to steer to a wide variecy cf trajectory profiles, 
such as velocity vs altitude profile, acceleration vs altitude 
profile, etc. This option is particularly useful for reentry 
trajectory shaping. 
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Generalized Acceleration Steering. - Under this option the 
steering angles are computed by solving a set of user specified 
equations. The dependent and independent variables in these 
equations must be selected from the dictionary of variables al- 
ready computed in POST. The only restriction is that these equa- 
tions must be explicityly a function of some derivative computed 
in the inner loop of the program. As a consequence, this option 
cannot be used to solve equations that are functions of integrals 
of the equations of motion. For example, this option cannot be 

used to maintain constant altitude by zeroing h. This is because 
the time derivative of altitude is computed from velocity, and 
velocity is computed from the integral of acceleration. The lin- 
ear feedback model should be used to solve problems involving in- 
tegrals of the equations of motion. 

In more precise terms, the steering variables are determined 
from the iterative solution of the problem: 

For each instant of time, determine the values 
of the steering variables, that satisfy the 
steering equations, 

e(9) *2.(1) " ^ - 0, (V-66 ) 

where ^ is a n-component vector of dependent variables, ^ is 

the desired value of these variables and £ the error in dependent 
variables. 

Typical applications of this option are given as 

1) Control normal acceleration to lg and axial acceleration 
to 3g by calculating the angle of attack and throttle 
setting that solves the equations 

A xB n) “ 96 *6 - 0 (V-67) 

A zB (a » n) - 32,2 " 0 

2) Obtain level unaccelerated flight by solving the equa- 
tions 

(a » n) “ 0 (V-68) 

Y a (<*, n) - 0 . 

Level unaccelerated flight is implicitly acieved in example 2 be- 

• , 
cause ■ 0 implies constant velocity, and m 0 implies con- 
stant altitude (that is, if y a ■ 0 when this option is initiated). 
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VI. TRAJECTORY SIMJIATlON 


The following sections present tie equations used in the 
trajectory simulation subroutines. These equations summarize the 
principal computations performed by tne program, and motivate many 
of the program input procedures. 


Events/Phases 

Simulation data are input according to phase, where the 
phases are defined by a user-specified sequence of events. The 
simulation equations are then solved sequentially by phase. 
Therefore, the user is required to input a sequence of trajectory 
segments that define the problem being simulated from beginning 
to end. These trajectory segments, or phases, are defined by 
two events — a beginning event and an ending event. An event is 
an interruption of the trajectory simulation that occurs when a 
user-specified variable reaches a user-specified value. An event 
must be created whenever the user wishes to change any input data 
for the problem or to cause any change in the method of simulating 
the problem. For example, the sequence of events for a typical 
ascent problem could result in a simulation setup similar to that 
shown in figure VI-1. 



Initiate pitch rate 4 20 aec after staging 

Initiate yaw rate 1 100 sec after event 7 

Orbit injection at inertial velocity of 25 568.0 fps 


Figure VI-1.- Event Sequence Setup 
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The event numbers for a given problem must be specified as 
real numbers by the user in monotonic increasing order. These 
event numbers are then used by the program to determine the order 
in which the events are to occur. The program requires that each 
problem have a minimum of two events — an initial event and a final 
event. Since a phase is initiated by the corresponding event, 
the event criterion for a given event specifies the conditions 
at the beginning of the corresponding phase. A problem is termi- 
nated by specifying the last event that is to occur. The problem 
can also be terminated in a psuedo-abort mode by specifying the 
maximum trajectory time, maximum altitude, or minimum altitude. 

Although event numbers must be monotonic increasing, they 
need not be consecutive. This allows the user co easily add or 
delete events from an input deck. 

Four types of events have been defined to provide flexibil- 
ity in setting up a given problem: 

1) Primary events - These describe the main sequential 
events of the trajectory being simulated. These events 
must occur, and must occur in ascending order according 
to the event number. Most problems will usually be sim- 
ulated by a series of primary events; 

2) Secondary events - These are events that may or may not 
occur during the specified trajectory segment. Secon- 
dary events must occur in ascending order during the 
interval bounded by the primary events. The occurrence 
of a primary event will nullify the secondary events 
associated with the previous primary event it they have 
not already occurred; 

3) Roving primary events - These events can occur any time 
after the occurrence of all primary events with smaller 
event numbers. They can be used to interrupt the tra- 
jectory on the specified criterion regardless of the 
state of the trajectory or vehicle. 

A) Repeating roving events - These events are the same as 
primary roving events except the interrupt values are 
input differently. There are two options for criteria 
value specifications. Option 1: Input the initial value, 

the increment, and the number of times the event is to 
be repeated. Option 2: Input an array of event cri- 

teria values. 

The program monitors as man> as ten events at a time, depend- 
ing on the types of events to determine which event is to occur 
next. This gives the user a powerful tool for simulating complex 
problems. 
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Multiple events are 


monitored in the following sequence: 


1) The next primary event is monitored; 

2) As many as nine primary roving events are then monitored, 
provided there are no secondary events. A roving pri- 
mary event is added to the list of those being monitored 
as soon as the primary event immediately preceding that 
roving event has occurred; 

3) Next, as many as nine secondary events are monitored, 
provided there are no primary roving events. (Note that 
caution must be exercised when using secondary events 
because of their nature. Since as many as nine sec- 
ondary events are monitored at a time, any one of those 
nine will occur as soon as its criterion has been met. 
Because they are secondary events, the event that occurs 
will cancel all secondary events with smaller event 
numbers .) ; 

4) Finally, a total of nine primary roving and secondary 
events: are monitored. 


Since the program can only monitor nine events (in addition 
to the next primary event), the sum of the primary roving events 
and the secondary events must be less than or equal to nine or a 
error will * jsult. 


The time-i-o-go model (TG0M) determines when the events 
occur during the trajectory simulation. Basically, TG0M checks 
T. valuU of the criteria being -nitored at “ 

step. If none of the criteria values has bracketed the desired 
cutoff value, then another integration step is taken. It a 
criterion variable is bracketed with the input step size, then 
TG0M computes a new stepsize equal to the predicted time-to-go. 

The predicted time-to-go for each event is computed from the 
equation 

At* - - y(t) At/ (y(t + At) - y(t)) 

where y(t) is the difference between the actual and the desired 
value of the event criterion. If more than one event is bracketed, 
then the minimum predicted time-to-go is used as the integration 
stepsize. This process is repeated until the criterion value 
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within the specified tolerance of the desired value. If the 
desired condition cannot be achieved in 20 iterations, an error 
message is printed and the program stops. Generally this situa 
tion is caused by an input error. The fundamental features of 
the time-to-go logic are shown in figure VI-2. 


Figure VI-2.- Illustration of Time-to-Go Logic 


Translational Equations 

The translational equations of motion are solved in the 
planet-centered inertial coordinate system. These equations are 


r = V 
L l -I 


i, ■ I ""' 1 [*r B ** A8 ] + £r 


where A is the thrust acceleration in the body frame, A 

-AB 

is the aerodynamic acceleration in the body frame, and Gj is 
the gravitational acceleration in the ECI frame. 


Initialization .- There are five options for initializing 
the velocity vector and two options for initializing the position 
vector. These options are described below. 

Inertial position components j Xj # , z \)'~ ^ e inertial 

position components can be input directly since no transformation 
is required. 

Earth-relative position or 6, ^ or h or rj.- this 

option the equations vary and the sequence of calculation varies 
according to the choice of input. However, the basic equations 
used are: 


9 I " 9 + % ( c - c o) 


(j> = tan -1 | k 2 tan <*> 


h + R s l* e 


if 6 is input. 


if <(> is input, 
g 


if h is input. 
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cos 4> cos 0 T 
c i 


r = i' cos <fi sin 6 
-II cl 
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Inertial velocity components |V^, ^yi* These variables 

can be input directly. 

Inertial local horizontal jV^, A^j.- * nert * a ^ com ~ 

ponents in the horizontal frame are first transformed to the 
geographic frame as 

COS Yj cos A zil 

— IG “ V I COS y I 8ln A ZI 
-sin y. 
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and then transformed to the ECI frame by 





„ Instantaneous velocity additions - A t the beginning of each 
” !. an instantaneous velocity can be added in the direction of 

i!lnuJ h H? 8t y® Ct0r * The "“gnitude of the velocity addition can be 
in P ut directiy, or can be calculated from the rocket equation 
trom the amount of propellant consumed 


AV " *0 * n 


r w i 

W G 

W - W 

LG W PC J 


If .. 4V is input, then the amount of propellant required to 
achieve this change in velocity is given by 

W pc - “ c [i - «P(- Av/e 0 I sp | 

The inertial velocity after the impulse is then given by 

*cos i cos i 
y P 

Yj - V~ + AV (IB] -1 sin i 


L cos i sin i S 
v p J 

° ptlon is 8*nerally used to simulate short burns for 
orbitai maneuvers. The direction of the impulse is controlled 
via the attitude of the vehicle and the engine gimbal angles. 


Static Trim 

an*lM e o r?ki%f rl V^ l0n 18 U8 * d to calcul «te the engine gimbal 
Jifoh h flaP deflectlon angles required to balance the 
pitch and yaw moments cuased by the thrust and aerodynamic forces. 

The static moment equation, 

M - 0, 

j! ® e T a11 ? nonlin «ar, and thus an iterative algorithm is used to 
th « solutions. This algorithm Is a suocm^ 



and the moment for the trimming engines is 


T M2P ■ c °‘ S E V co ‘ 

ieN T 1 



Az 

P ? 


sin i 

P 



C08 



The static moment equation in the pitch plane is then given by 
**P “ *MP + T M1P + T M2P* 

where A^p depends on the flap deflection angles and T^p 
depends on the engine gimbal angles. 

Moment equation in the yaw plane .- The aerodynamic moment and 
the thrust moment for the i*b engine are shown in figure 17. 
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The total aerodynamic moment in yaw is 

V - 1 S ( C n *RY + ®r ' C A ly R)' 

The total thrust moment in yaw for the nontrimming engines is 


M1Y 


L T / cos i cos i Ay + sin i Ax 

R i \ P i y i P i y i P ±) 


i^N_. 
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and the moment lor the trimming engines is 


II2Y 


■ cos i 


L T cos i Ay 
K P, e 


ie»„ 


+ sin i V T_ Ax 

y \ p 4 

ieN_ 
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The static moment equation in the yaw plane is then given by 

+ T M1V + W 

depends on the flap deflections and ^J2Y depends on 
the engine gimbal angles. 
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Integration Methods 

The number of integrals computed during any particular phase 
is determined from the options requested by the user. As a 
minimum, the translational equations of motion ere integrated to 
give the position and velocity of the center of mass of the 
vehicle. The user may also select additional variables to be in- 
tegrated. The only restriction is that no more than 30 integrals 
can be computed per phase. 






POST contains 3 general purpose Integration methods and sev- 
eral special purpose orbital propagation methods . These methods 
are summarized below. 


Runge-Kutta Methods .- POST contains two Runge-Kutta integra- 
tion methos: (1) the standard 4th order method, and (2) a 8th order 

method. The calculations for these integration methods are based 
upon the formula 


my n* £ b l k i’ 

where 


\ - hf/: 


x + 


V 


Cjh, 


s 

+ L 

i-l 


‘u k j) • 


i“l»2 ,"',s. 


These formula are 


represented by the array 


C 1 

*11 

a 12 

* * * a ls 

C 2 

• 

*21 
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a 22 
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a 2s 
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• 

*sl 

• 

• 

®s2 

• 

• 

■ • • & 

88 


b. 

b. 

• • • b 


1 

2 
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The coefficients for the two methods are given in Tables VI-1 and 
Vl-2, respectively 


Table 1V-1 

Runge-Kutta 4th order (Kutta) 





Runge-Kutua 8th order (Shanks (1966)) 


0.000 
0.111 
0.167 
0.250 
0.100 
0 . 167 
0.500 
0.666 
0 . 333 
0.833 
0.833 
1.000 


0 . 11 U 

0.0042 

0.0069 

0.0580 

0.0340 

- 0.5838 

- 0.1235 

3.6265 

0.9043 

0. 1250 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

0.1875 
0.0780 
0.0000 
0.0000 
0 . 0000 
0.0000 
0.0000 

- 0.0360 
0.0041 
2.1111 
- 0.1317 
' - 10.6667 
- 2.6296 

0.1286 

3.4722 

0.5144 

19.2901 

4.2438 

- 4.5000 

0.0000 

26.0000 

5.6667 

0.4074 

0.7469 

- 0.3642 

- 0.0833 

0.5000 

1.000 



0.8043 

0 . 0000 

0 . 0000 

2.6296 

- 4.2438 

6.1667 

0.6358 

0.0000 

0.0000 

0.1000 

— 

- 1.9411 

0.0000 

0.0000 

6 . 937 ? 

11.0095 

- 14.9268 

0.0854 

- 0. 1646 

- 0.4390 

- 0.2927 

0.7317 

0.0488 

0.0000 

0.0000 

0 0.0000 

0.0000 

0.2571 

0.3238 

0.0321 

0.0321 

0.0429 

0.2143 0.0488 


Krogh Variable-Step Variable-Order Integrator 

Tne variable-step length variable-order predictor-corrector 
developed by F, T* Krogh of the Jet Propulsion Laboratories pro- 
cedure represents the state-of-the-art in numerical solution of 
systems of ordinary differential equations. It includes all of 
the following facilities: 

1) A core integrator for advancing the solution from one 
uniform step to the next consisting of variable order 
Adams predictor-corrector equations requiring the stor- 
age of a difference table for only the highest ordered 
derivatives; 

2) A method to scare integration with first-order equations 
and increase the order to as high a level as numerical 
stability permits; 

5) Algorithms for changing the step size and updating ac- 
cordingly the difference tables of the highest-order 
derivatives including appropriate smoothing to prevent 
numerical instability; 

4) Algorithms for deciding when and by how much to change 
the step size baaed upon the accuracy requested by the 
user; 

3) Tests for numerical stability of the predictor-corrector 
order and step size tentatively chosen in the context of 
the current differential equation set; 
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Test for accuracy requests that are so stringent that 
round-off error prevents their satisfaction; 

7) Algorithms for tihe automatic selection of the core inte- 
grator to fit the characteristics of the set of differ- 
ential equations at hand. 

8) An interpolation algorithm for obtaining dependent vari- 
able values to the user-specified accuracy at values of 
the independent variable different from normal integra- 
tion steps. 

The Krogh integrator was developed to meet the conflicting 
goals of (1) reliability, (2) efficiency (in the sense of mini- 
mizing the number of derivative evaluations to obtain a given 
accuracy), (3) flexibility, (4) low integration overhead, and (5) 
email storage requirements. The goals are listed in the order in 
which they are emphasized in the procedures. The package lias no 
equal in reliability. All the user need provide the integrator 
is his accuracy requirement and a tentative step size that is 
less than the integration interval. The integrator then uses 
all eight of its algorithms to provide a solution requiring as 
few derivative evaluations as possible while remaining within 
accuracy tolerances. 


This documentation restricts itself to algorithms (1), (2), 
and (3). Algorithm (3) is thoroughly described in Reference 9. 
Algorithms (4) through (7) are not amenable to rigorous mathe- 
matical analysis, but rather have evolved from extensive numer- 
ical experimentation by Krogh and his associates at JPL. They 
are adequately described in the source code and references. Fur- 
ther, to simplify exposition, equations will be given for the 

single d ^ order differential equation 

Y (d) - f(t, Y, Y 1 , •••, Y d_1 ) [1] 

The following notation is used: 


l - independent variable 
h - integration step size 

y(t) - corrected (final estimated) value of Y(t) 
t - value of t at current integration step 


t 


n+k 


t + kh 
n 


[ 2 ] 

13] 

14] 
[5] 

16] 


v v v ••• 


y(t n ), Y(t n ), y 1 ^), *•* 


[7] 


predicted (Initial estimated) value of Y(t ) 

(d-l)l 


f ■ f/t » y » y*» ••• • y ^ ~ ') 
n | n 7 n n 7 n ' 


[3] 

[9] 


V! -14 


v 1 * 


v x-1 f - v 1 " 1 f . 

n n-1 


for i-0 

for i*l,2, ••• 


t f f £ (v V •••« p » d ' 1) ) £or i '° 

P D ( 4-1 <-1 


[ 10 ] 


[ 11 ] 


V i-1 f - v 1 ” 1 f 
P n n-1 


for i=l,2, 


The predictor and corrector equations of the core integrator 
are obtained by integrating Newton's interpolating polynomials 
for the most recent q and q+1 values, respectively, of the highest 
order derivative, f. For the predictor, f is approximated as the 
(q-l)th degree polynomial 


q-1 


£(X n +Sh> ' 


r 12] 


i-0 

Successively integrating approximation [12] from s-o to e-l yields 
the predictor equations. 


P w-j) -y h k , (d -i*> + hJ 

n+1 rr n 

k-0 

q-1 

^ Y A j V i f n for j-1,2, • • • ,d 
i-0 

[13] 

wh,r ‘ 

for 

i*i, * * • ,q» 

[14] 

d ij(*) ~5 0 (t) dr 

for 

i-i,***»q. 
j“l» * * * ,d 

[15] 

and Y^j * o^d) 

for 

i-i. ,#, .q 

J-l* , * , .d 

[16] 


For the corrector, f is approximated as the q th degree poly- 
nomial 

f c • f** +<j /a] V*f <y f (x + sh) . 

n+s n+s qO w p n 


[17] 
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Successively integrating approximation [17] from s "0 to s*l gen- 
erates the corrector equations 


y 


(d-j) 

n+1 


P 


(d-j) 

n+1 





for j«l, •••,<!. 


[ 18 ] 


Notice that in the equations [13] and [13] the coefficients Y^j 

are independent of the order, q. Hence, Adam f s method has the 
practical advantage that its order can be increased simply by 
adding more terms to the predictor and corrector equations with- 
out revising their coefficients. Systems of differential equa- 
tions are handled by applying equations [13] and [18] to each 
component of the system. Different components may have differ- 
ent values of d and q. 

The starting algorithm requires no special techniques such 
as Runge-Kutta equations or Taylor-series expansions. Instead, 
q i 3 simply taken as 1 in the normal predictor and corrector equa- 
tions. There is, nonetheless, i starting phase during which the 
highest-order derivative evaluation following the calculation of 
the corrected values of the integrals is skipped. During this 
starting phase, the values of q and h increase quite rapidly. 

When numerical stability problems appear during this initial phase, 
the ordinary corrector equation, [18], is replaced by 


ur j) ■ p nir j> + hJ - *ij ^ £ n+l iot j- 1 - 2 -— 

Here a is chosen on the close > interval from \ to 1 to give opti- 
mal stability characteristics. 


[19] 


Algorithms [3] for changing the step size are only equipped 
to double or halve h. Hence, algorithm [4] for selecting the 
step size is restricted in the same way. Both procedures are, 
however, very efficient so that adjustment to the optimal step 
length occurs as quickly as reliability permits. The two algor- 
ithms are described in detail in Reference 11. 


Like the predictor and corrector equations, the interpolation 
equations of algorithm [8] are also generated by successively in- 
tegrating a degree interpolating polynomial for the highest 
order derivative. Indeed 





5 : f(x + sh) 
n 


[ 20 ] 
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[21] 


(d-1) . (d-1) + v P U 

4 * a ^ n 4 —* 


7 ‘ C. * Y<d ' 1) <*„ + •» 


n+s 


4-0 


where 


il 


-j; <■ r)“ 


Evaluating the differences* V ^ n + g * U3 ^ n 8 definition [20] and 
substituting the results into definition [21] yields 

q 


y - y^ + V 

*n+s y n Lu 


i-0 


4-0 


V 1 f 


n+1 


[ 22 ] 


[23] 


Generalizing this procedure to derivatives of lower order produces 

j-i q-j 

t (d-j+k) + h j 

rr n 

i-0 

for j-1, ‘‘‘.d 


■£ 3) •E^. Wni E 

k-0 fe * 


[24] 


where 




[25] 


and 


q-j 


r i ^ ( * ) ^ P 41 T (i-i+l) (j-1) 

i-0 


H • 


[26] 


Laplace' s method.- Laplace's method is an iterative technique 
for propagating thVposition and velocity (in planet-centered 
coordinates) of a nonthruating vehicle in vacuum during flight 
over a spherical planet. The technique is based on the analytical 
solution of the two-body equations, and yields the inertial state 
at time t + At as a function of the state at time t. j 

The equations used in Laplace's method are: 

r (t + At) » f iTj (t) + g V_x ( fc ) 

V x (t + At) - f r x (t) + g Vj (t) , 

where the scalar coefficients f, g, f, an< * S ar ® gi ven 1° 
Table VI -3- 
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Table VI- 3 - f and g Series 


Coefficient 

* 

Elliptical 

* 

Hyperbolic 

f 

. ^ a 2 AE 

1 - 2 — 8in z — 
r 2 

n 

1 + 2 — sinh 2 

r 2 

g 

At - ~ (AE - sin AE) 

At - — (AE - sinh AE ) 

• 

- /pa sin AE 

-/ua sinh AE 

f 

r r . i 

n n+1 

r n r n+l 

• 

g 

a 2 AE 

1 - 2 - sin^ t— 

1 + 2 — sinh 2 ~ 

r n+l 2 

r n+l 

*Note that 

r^ - li(t + j At) | , and n * ~ . 


As indicated in the above table, the f and g coefficients 
are computed analytically from the change in eccentric anomaly 
during the time period At. This is in contrast to the standard 
Laplacian method, where these coefficients are infinite series 
in At. The change in eccentric anomaly is calculated by solving 
a special version of Kepler's equation via the Nevton-Raphson 
algorithm 


AE 


a+1 


F^AE) 


1 , 2 , 


10 . 
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This algorithm converges rapidly from the initial estimate 

AE Q - n At. (VI- 29) 

The fora of Kepler’s equation that is computationally ef- 
ficient for this application is 


sin AE + 2 


— sin 2 ~ + AE - n At , 
ni a for e < 1 


r • v 

sinh AE - 2 sinh 2 AE + AE - n At, 

^ a for e > 1. 


The derivative dF/dE used in the Newton-Raphson algorithm 
is given by 


F'(AE) 


([H 
1 1 ?-] 


r • v 

cos AE + — — sin AE f 1, for e < 1 

i /y a 


r • v 

cosh AE - — — sinh AE - 1, for e > 1. 

a 


Encke ' s method .- The Encke method u3ed in this program is 
modified from the usual Encke technique in that it rectifies the 
reference conic at every integration step and does not use the 
standard Q-series expansion in calculating the gravitational in- 
crement . 

The Encke method should be used for orbital problems where 
small perturbing accelerations, such as the oblateness of the 
planet, atmospheric drag, or solar electric propulsion, must be 
included in the simulation. Numerical results indicate that, for 
problems involving small perturbations from Keplerian motion, 
Encke 1 s method is approximately four times faster than Cowell 
methods, which integrate the total acceleration. 

The Encke method determines the total motion by summing the 
motion due to the two-body equations and the motion due to the 
perturbations to the two-body equations. The position and veloc- 
ity in the inertial planet-centered system at time t + At is 
given by 



4 


i 

. i 

j < 

j ■ 


i 

I 

i 

\ 
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iTj (t + At) « £2 (t + At) + Ajr (t + At) 
Vj (t + At) - V 2 (t + At) + AV (t + At) , 


(VI- 32) 


where (£ 2 , V? ) denotes the Keplerian motion computed by Laplace's 
equations; that is. 


£2 (t + At) - f £2 (t) + g V2 (O 

V 2 (t + At) - f £2 (t) + g V2<t), 


} 


(VI-33) 


and (A£, AV) denotes the numerical solution of the differential 
equations 

Ar - AV 

« • [[“!■* [in, + *„] + £,] - 12 (a + Sr) 

Ar(t) - AV (t) » 0 t 

where £2 (.£2 + ArJ is th® two-body acceleration at r? 4- Ar 
The Runge-Kutta or Adams -Bashforth method can then be used to 
the above equations. 

Integration step models.- The integration step t At, is gen- 
erally specified in terms of an increment in time. However , this 
option enables the user to specify the integration step in terms 
of an input increment in true anomaly. This option is useful in 
orbital problems where the geometry is easily expressed as a func- 
tion of true anomaly. 


(VI-34) 


solve 


The following equations are then used to calculate At as a 
function of A0. 


02 * 9 } + A0 

£ 2*2 tan 


1 |/_l-e\ 2 

I \ 1+e J 


tan 


r 4 . - a(l - e cos E 2 ) 
A£ * E2 - Ei 




At • 

r 


(AE - sin AE) + r2 Ii sin A0. 
H 
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In these equations, subscripts 1 and 2 denote current and future 
values, raspsctively. 


-2i 
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Launch Options 


There are two specific options for simulating particular 
launch conditions; (1) hold-down for vertical takeoff , and (2) 
hold-down for horizontal takeoff* These specialized options are 
required to simulate certain physical constraints that are not 
modeled in the equations of motion. 

Holddown for vertical takeoff*- This option is used to simulate 
vertical (rocket type) takeoff. When using this option, the rel 
ative position and velocity remains constant while the inertial 
position changes by the Earth rotation. The inertial velocity 
magnitude remains constant while its direction changes. This 
model simulates physical constraints that hold the rocket on the 
launch pad until the rocket is released* The equations used to 
calculate the accelerations that produce this motion are 


* 

4 


Aj - £ x V 


(vi-36 ) 


Holddown for horizontal takeoff .- This option is used to simulate 
horizontal (aircraft type) takeoff. When using this option the 
vehicle accelerates in the local horizontal plane according to 
the forces described by the user input. The vertical component 
of acceleration is internally computed to produce the proper hor- 
izontal motion. The equations used are: 


- A -*c - [“] 


^SI' 


(VI-37 ) 



VII. AUXILIARY CALCULATIONS 


In addition to computing the basic variables, POST also com- 
putes numerous auxiliary variables that are related to: (1) 

conic parameters, (2) range calculations, (3) tracking data, (4) 
analytic impact calculations, (5) velocity losses, and (6) veloc- 
ity margins. The equations used to calculate these variables are 
presented below. 

Conic Calculations 

The following Keplerian conic variables are computed. 

v T 2 

C I U 

& energy, -z — 

L I 

a seraimajor axis, -u/ 2& 

h angular momentum, |xj x 

p semilatus rectum, h 2 /u 

e eccentricity, /|1 - p/a| 

AV velocity required to circularise orbit, /AV • AV, where 

^ -h/h 

Sri ' V r I 

S, " Sh " WK * Sri 1 

y,. - (u/t,)'* ^ 

AV * V - V T 
— -c — I 

i inclination, cos" 1 (h^/hj 

ft longitude of ascending node, cos -1 (x^ • u^)> where 
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V 



Ufl - Zj x h/llj x h| 


argument of vehicle, p * cos” 1 | i 


time since perigee, ^ M 


time to perigee, P - T g p 

latitude of perigee, tan" 1 [u$/ + uf} , where 

u - cos(o))u^ + sinCo))^ x 

longitude of perigee, tan” 1 (U 2 /U 1 ) 


altitude of perigee, r - R (<J> ) 
r ° p s p 

altitude of aoogee, r - R ($ ) 

a s p 



Z 1 + 

_e\ 

a 

U- 

e ) 


velocity at apogee, J j 
hyperbolic excess velocity, / 2 6> 

maximum true anomaly for hyperbolic orbit, cos -1 (-1/e) 
declination of outgoing asymptote, sin" 1 [u r (3)], where 


u r ’ CO9(0 max " e) Sri + 8ln (6 max _ 6) % 


R. right ascension of outgoing asymptote, tan* 

A 


C'^ 
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e 


E 


true anomaly, cos -1 



eccentric anomaly. 


2 tan" 



- e 
+ e 



) 


M mean anomaly, E - e sin E 


0 ) 

r 

P 

r 

a 

P 


argument of perigee, p - 6 
perigee radius, a(l - e) 
apogee radius, a(l + e) 


period. 



Range Calculations 


The progam provides for various types of range calculations. 
The equations for these calculations are given below. 


Dot product downrange .- The relative range angle, measured 
from the vehicle's initial position to its current position, is 
given by 



(VII-1) 


where ii is a unit vector along the initial position vector 
so 


in Earth-centered rotating coordinates and u is a unit vector 


along the current position vector in Earth-centered rotating co- 
ordinates. The range over an oblate spheroid is calculated from 
the average radius to the surface, and is given by 



(VII-2) 


Crossrange and downrange via orbital plane reference .- Re- 
ferring tc figure VII-1, identify the vehicle's position at time 
* * 
t by 0, and at a later time t by P. At time t , the ve- 

* * 

hide has a latitude of $ , a longitude of 6 , and a velocity 


! 


3 i 

ii 

j 

j i 


I 


■j 

J 
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heading of X . At time t the vehicle is at latitude ♦ and 
longitude- 6 . The downrange angle (y) and the crossrange angle 
(v) shown in the illustration are measured along, and normal to, 
the great circle through 0, and are inclined to the meridian by 
X*. From analytical geometry, v and y can be expressed as 


sin v - -sin X* sin cos <J> c cos 6' - cos X cos <|> c sin 0 


* * 

+ sin X cos $ sin 


a ^ * 

sin y - (-cos X sin $ cos <J> cos 9' + sin X cos 4> c sin 0 


* * , \ i 

+ cos X cos <j> sin <J> c j/cos v 


(VII-3) 


( it & v . 

cos ^ cos <(> cos 0 + sin <£ sin 4> \/cos v, 


where 9* and X can be defined in either of two ways: 


1) The great circle to which v and y are referenced is 
fixed and rotating with the Earth. Then 


X* * Earth's relative heading - sin -1 


u^ + v 2 


R ' R 


9" «" e - 0 ; 
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2) The great circle to which y and v are referenced is 
inertially fixed, having the Earth rotating below it. 
"hen 


X* - inertial heading - sin -1 


u 2 + v 2 


(VII-5) 


e'«e-0*+fl (t - t*). 
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V 


t 

\ 




Knowing v and y, and crossrange C R and downrange D.. 


distances are 


C„ - R v 

R ave 


D„ « R y, 
R ave 


where R Is the average Earth radius between the Initial and 
ave 

final points. 


t 


i 
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Auxiliary Position and Velocity Calculations 


The solution from the translational equations is then used 
to calculate numerous output variables. The key variables directly 
computed from (Xj, y^» Zj) and (V X j, V ZI ) sumarized below. 

Tj “ geocentric radius 

' ( £ i ' Li t \ 

Vj ■ magnitude of the inertial velocity 

■ (2i ' h)' > 


relative velocity 


V — Q x r 
-i -i 


* atmospheric relative velocity 


V + V 
-R -W 


magnitude of the relative velocity 

i \t .1? ^ 


/ v ^ 

2 r - 2 r)' 


* magnitude of the atmospheric relative velocity 

- (v. • 

\ — A —Ai 

u^j ■ unit vector along radius vector 
’ ^l/ r i 

Uyj ■ unit vector along inertial velocity vector 

■ V v i 


Yj “ inertial flight path angle 


» sin" 


rll 
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r 


f 

[ 

r 


l 



'R 


relative flight path angle 
8 In" 1 


— IG 


^RG 


v 

-AG 




*ZR 


ZA 




atmospheric relative flight path angle 

sltr ‘ [2 h • Sva] 

Inertial velocity in the G-frame 


[IG] Vj 


relative velocity in the G-frame 

' i6 >2r 

atmospheric relative velocity in the G-frame 
[IG] V A 


inertial azimuth 

[ v yg/\g] 
relative azimuth 

tan 1 [ v ryg/ v rxg] 

atmospheric relative azimuth 

tan 1 [ v ayg/ v axg] 

geocentric latitude 
•in-> 

inertial longitude 


i 

< * 
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relative longitude 


9 I - «p I* - c o) 

sensed acceleration In the B-frame 
^TB + -AB 

magnitude of the sensed acceleration 



sensed acceleration in the ECI-frame 


Auxilary Attitude Calculations 

The attitude angles that are not used to generate the steer- 
ing commands are computed for output in the auxiliary calculation 
subroutine. These equations are summarized below: 

1) Aerodynamic angles : 


tan 1 ( v azbAiXb) 

ta "-‘ (Vb//' 

tan -1 ( 


2 + 
AXB 


V 2 \ 
AZBl 


) (VII- 8 ) 


GB 23 + sin 6 sin^ 


GB 2 2 c <>8 - GB 21 sin * 2 * CO * Y A I 


2) Inertial Euler angles: 


■ tan -1 (LB 23 /tB 22 ), 
*! “ -sin ’ 1 (LB 21 ), 

■ tan_1 (LB 31 /LB n ); 
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3} Relative Euler angles: 
V R - tan ' 1 (GB 12 yfcBn), 

0 R - -sin-^GBu), 

♦ R ■ tai |GB 23 yiGB 33 j . 


) (VII-10) 
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Tb'' relationship between the body rates and the attitude 
angles are: 

1) Aerodynamic angles : 


w I a l + + ( sin °0 B 


w . - • 2 + d ? o + a 


a > 2 1 a 3 + d 30 - (cosa) 6 


2) Inertial Euler angles: 


t>j. cos ^ cos 0 ^ - sin 0 ^ 


w y I " ®I ‘ *1 8in *1 


► cos ip sin 0 + \p cos 0, 

I 1 I L 1 


3) Relative Euler angles: 


>R - *R 8l " ®R 


w - [GB ] — 

y r 


+ e R cos + ip R sin <t> R cos Q R 


tan <t> 


i R cos $ R cos 0 R - e R sin <J> R 
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Tracking Data 


POST computes tracking information for as many as ten 
tracking stations per phase. The tracking stations are located 
on a reference ellipsoid and are specified in terms of their lati- 
tude, longitude, and altitude above the ellipsoid. These variables 
are illustrated in figure VII-2. 



Figure VII-2.- Radar Tracking Schematic 
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The position components of the tracker in the Earth-relative 
frame are given by 

cos * T cos e T 
cos * T sin e T 
sin ♦_ 

where is the altitude of the tracker, $ the latitude of 

the tracker, and 8 T the longitude of the tracker. 

The slant range vector in the ECI frame is given by 

ISR “ -I " IIPl_1 — TR* (VI 1-15) 

and the slant range is then computed as 

r SR " >/- SR * ISR* (VII-16) 

The elevation angle can then be computed as 

t t • ,ln ' 1 (are • r -s R / r SR) ■ (VXI-17) 

where 

%R " ^Tr/ I IIPl " 1 -trI • (VII-18) 

The slant range vector, transformed to the geographic frame, 
is 

^SRG “ tIGl ^SR’ (VII-19) 

and thus the tracker's azimuth is given by 

*2t ■ (^srg/ *srg) * (VII-20) 

The look angles are calculated from the slant range vector 
transformed to the body frame; l.e., 

IfiRB " (IBl (VI1-21) 
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Using the coaponents of £ggg» the cone angle Is then given by 

*T - co.-> (* siJ ,/r SR ), (VII-22) 

and the clock angle Is given by 

°c ■ “”■* ( y S»/*Sl»)' <»«-») 

Space losses are calculated for the tracking stations as 
follows: 

SLj - 36.56 + 20 Log l0 (R slm • FR^ 

SL 2 - 36.56 + 20 Logjo ( R SLM * FR 2 ) V (VII- 24) 

SL 3 - 36.56 + 20 Log 10 (R gLM • FR 3 ) , 


where 


FRi 


420*0 (command frequency) 


^2 « 2287.5 (telemetry frequency) 

fR 3 * 5765.0 (tracking frequency) 


* slant range distance in statute miles, 
M 


Analytic Impact Calculations 

The analytic impact calculations predict the geodetic lati- 
tude, longitude, and time of flight at impact for a vehicle with 
a given position and velocity to Its intersection with the sur- 
face of the oblate planet. These calculations assume Kepler ian 
motion and are no" corrected for drag effects. 

The basic problem in determining an impact point from a 
specified position and velocity jXjo* -10) * 8 in ca l cu l atin 8 

the impact eccentric anomaly. This angle is determined by 
iteratively solving the equation 
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r 


r i< E) • M*c) + V 

where h is the desirel impact altitude above the oblate 

ip 

planet and the position vector is giveu by 
r r (E) - f(E) r IQ + g(E) V IQ 

f (E) - (cos (E - E q ) -e cos E Q ) / (l - e cos E Q ) 


-JY < 


sin E - E rt | -e sin E + e sin E, 


Once the impact eccentric anomaly, E^ , is determined, then 
the time, latitude, and longitude of impact are calculated as 

t « t„ + I — /E. - E. -e sin E. + e sin E-.\ 

ip 0 J y V ip 0 ip 0/ 

a 

♦, lp ' + y ip) 

0. ■ tan- 1 ( 0 t. 

\ x ip ) p lp 
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Velocity Losses 


There ere two velocity loss options: (1) inertial velocity 

losses , and (2) atmospheric relative velocity losses. The iner- 
tial losses are used £or orbital problems while the atmospheric 
relative losses are used for atmospheric flight problems. 

The total change in the magnitude of the velocity is given 
by: 


V f " 


AV - AV* - AV__, - AV 


TV 


AERO - 


AV - 

g 


AV 


atm 
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where 


iV. dt 

4 V -f: s) Jt 

•K 


AERO 


F. • U 
-A — 


AV 


-f 


g sin y dt 


AV 


atm 




dt 


ideal velocity 
thrust vector loss 
aerodynamic loss 
gravity loss 

atmospheric pressure loss 


Inertial losses are computed when IJ • [IBjVg/V^, and simu- 
larly atmospheric relative losses are computed when U * (IBJV^/V^. 


Velocity Margins 


The program computes the amount of velocity margin avail- 
able and the amount required , based on a percentage of the ideal 
velocity. 

The velocity margin is calculated as 


AV H ■ «0 I SF tn 


(V W G - W p) 


The excess velocity is then given by 


AV 


E 


a *m " AV % 
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Sun-Shadow Calculations 


The program computes several sun-shadow variables. These 
variables are used to calculate the sun-vehicle orientation angles, 
the sui-snadow conditions, and ."he position and velocity of the 
vehicle in the vernal equinox system. These auxilary variables 
are based upon the Greenwich hour angle (GHA), the Greenwich 
hour angle of the Sun (GHAS), the declination of the Sun 6 

and the time of reference past midnight (TRPM) . The Greenwich 
hour angle, the right ascension and declination of the Sun remain 
constant during the simulation. 


The vehicle position and velocity vectors in the vernal equi- 
nox system are given by 


•Ive 
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The Sun unit vector in the (VE) system is 

F*. c °.l 



8 


where 


(VII-33) 


- GHA + GHAS. 

The Sun unit vector in the ECI system is calculated as 
- [IV] 'u . 


The cone and clock angles of the Sun vector in the body system 
are given as 7 


6 

cone 


cos 


-1 




*clock " t4n ( y eB^*eB) 


i 
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where u fi is a unit vector in the direction of the axis, and 
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The program also computes a shadow function, which is used to de- 
termine when the vehicle is in or out of the Earth's shadow. This 
function is based upon a cylindrical shadow model and is given by 


where 
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a ■ V* • ± 

4-£l-(£l •«.!) 


If S<0, then the vehicle is in the shadow of the Earth. If s*.0 f 
then the vehicle is in the sunlight# The vehicle is entering the 

shadow if SlO, and exiting if S>0. 


4 


l 


Multiple Vehicle 

The program has the capability to simultaneously simulate 
the motion of two independent vehicles. One of the vehicles is 
active in the sense that it can be controlled using' propulsion 
and/or aerodynamic forces. The other vehicle is passive in the 
sense that ic cannot be controlled and is assumed to be out of the 
atmosphere and nonthrusting. As a result, the active vehicle is 
referred to as the pursuer, and the passive vehicle the target. 

The relative geometry between these two vehicles is defined in 
Figure VII-3. 

A large number of output variables are calculated for the 
target vehicle. These variables are computed using equations 
tha 1 . are identical to ':hose used for the active vehicle. A com- 
plete list of these aquations Is given in Volume II - Utilisation 
Manual. Only the key equations are in this section. 





VI 1-18 



The increment in position, velocity, and acceleration between 
the active and the target vehicle are given by 

A r - r, - r » Ar + / AV dt 

"° (V1I-37) 

AV » V ^ - v - AV 0 + / Aa dt 



Figure VII-3.- Relative Geometry between 
Active end Target Vehicles 
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VIII. TARGETING AND OPTIMIZATION 


POST uses an accelerated projected gradient algorithm (PGA) 
as the basic targeting/optimlzation technique. PGA is a combina- 
tion of Rosen's projection method for nonlinear programming (refs. 
3, 4, and 5) and Davidon's variable metric method for unconstrained 
optimization (ref. 6). The program also contains backup single- 
penalty function methods that use steepest descent, conjugate 
gradients, and/or the Davldon method. These standard gradient 
method are well documented in references 6 and 7 and are only 
brief 1> described in the following discussion. 

The projected gradient algorithm is an Iterative technique 
designed to solve a general class of nonlinear programming prob- 
lems. PGA employs cost-function and constraint gradient informa- 
tion to replace the multidimensional optimization problem by an 
equivalent sequence of one-dimensional searches. In this manner, 
it solves a difficult multidimensional problem by solving a se- 
quence of simpler problems. In general, at the initiation of the 
iteration sequence, PGA is primarily a constralnt-satlsficatlon 
algorithm. As the Iteration process proceeds, the emphasis 
changes from constraint satisfaction to cost-function reduction. 

The logic used to effect this changeover process will be dis- 
cussed below. 

Since numerous analytical developments of this technique are 
available (see refr 3, 4, and 5), this presentation will pri- 
marily emphasize the geometrical aspects of the algorithm. This 
geometric interpretation clearly motivates the equations and 
logic contained in PGA, and a basic understanding of these con- 
cepts is usually sufficient to enable the user to efficiently 
use the algorithm. 


Problem Formulation 

The projected gradient method solves the following nonlinear 
nrograonring problem: 

Determine the values of the Independent variables, u, that mini- 
mize the cost function (optimization variable) 

F(u), (VIII-1) 
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subject to the constraints (dependent variables) 

c(u) >0, (VIII-2) 

where u c R a ; c is a vector-valued function, i.e., £:R m -*• R n ; 

“ a 1 

and F is a scalar-valued function, i.e., F:R -*■ R . 

The algorithm is actually more vet’''*' lie than this simple 
formulation might indicate. In order to .maximize any particular 
function, say W(u), all that is required is to define 
F(u) ” -W (u) and determine the minimum of F(u). The equality 
constraint case is also contained within the above formulation 
since constraint equations of the form 

Cj(u) - 0 (VIII-3) 

are special cases of Sq (VIII-2) . 

In the trajectory optimization, the cost fu action and the 
constraints are not explicitly a function of the independent 
variables, but rather depend explicitly on the state variables 
ir , Vj, m, and Q. The explicity equations relating the state 

(dependent) variables to the independent variables are the in- 
tegrals 


ii ■ \ * fa ic 

f. 

m * ra + f u dt 

° J 

q - y*Q dt. 


dt 


(VIII-4) 
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denote* the above etate variables of the system being 


-n + 

simulated at the n event, and * 


and x ~ denote the value 


of x on the plus and ninus sides of that event, then 


-«*£. 4 


where u are the independent variables in phase n, and r 
—n 

represent the solution of the state differential equations over 
phase n. The values of the state variables on the positive side 
of event n are then 

x + ■ x + Ax (VI 

— n ~n n 

where Ax represents the discontinuity in state .(e«g«» velocity 

n 

impulse at the n event) • 

The cost function and the trajectory constraints are computed 
at the positive side of the specified events, and are therefore 
given by 


( U) - f(4 £ ) 


c(u) “ 


+ \"l 


'1 (s) 


[S K) J 

where v^ denotes the event at which the optimisation variable 
is specified and v^ denotes the events at which the dependant 

variables are specified. This generality enables the program to 
solve problems in which intermediate constraints are defined, as 
•fell as problems where the cost function is not specified at the 
final event. 
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The trajectory propagator, T , can represent either numer- 
ical Integration or analytical Keplerlan equations. 


Fundamental Concepts and Nomenclature 

To facilitate the discussion of the projected gradient algo- 
rithm, the following nomenclature and basic concepts will be In- 
troduced. 

A real k-dimensional Euclidean vector space is denoted by 

r\ and x denotes a column matrix whose elements are 
x ± i where” i - 1, 2 k. The vector inequality x >_ 0 im- 

plies x ± _> 0 for each i, and A" denotes the transpose of the 
real matrix A. 

The cost gradient is an m-vector of partial derivatives de- 
noted as VF or 3F/3u, and is defined as 

®>i ■ {^- < 

The gradient to the i th constraint is similarly represented. 

The Jacobian matrix of the constraint vector function with 

respect to the independent variable is a matrix whose i th row is 
the gradient vector Vc^. Thl® matrix is denoted as 

3c 

J(u) - -r- 0 


and contains n rows and a columns. Clearly, 


ij auj* 
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The J 1 * 1 constraint is said to be active at u if and only if 

a) c^u) < 0, (VII1-12) 

An active constraint is said to be unconstraining if and only if 

b) c <u) - 0 and r^ - [(SS') _1 £ 0. (VIII-13) 

Condition a) implies that the constraint is either vio- 
lated at u, while b) indicates that the negative of the cost 
function gradient "points" outside the feasible region. 

The sensitivity matrix is that matrix whose rows are the 
gradients to the active constraints, and is denoted by 

he 

S(u) " TZ* (VlII-14) 

where e is the n^-vector of active constraints. Equality con- 
straints are always active and thus are loaded into the upper 
elements of the e. Thus* £ is essentially the error vector 
for the active constraints. The error function is defined to be 

E(u) - e'e. (VII1-15) 

The sensitivity matrix, S, is obtained from the Jacobian 
matrix, J, simply by deleting those rows that correspond to 
inactive constraints. 

Corresponding to each constraint function c.(u) is a 
boundary hypersurface , B^, defined by 

B i “ |H sc i<“) m °J* (VIII-16) 

Clearly, is an m-1 dimensional nonlinear manifold. It can, 

however, be approximated at each point 0 in R* by an m-1 
dimensional linear manifold 

C^d)- iu:V Ci '(a) (u - d) + Cl (u) - oi. (VI 1 1.-17) 
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The feasible region for the inequality constraint is the 
half-space in the independent-variable space defined by the set 

L (VIII-18) 

while the complete feasible region for all of the constraints is 

n 

R - O R (VIII— 19) 

i-1 

The boundary of the complete feasible region must be 

n 

B(R) - U (B 4 n R) (VIII-20) 

i?l 

The intersection in the preceding equation is required to select 
from the unbounded boundary, of the feasible region of the 

l** 1 constraint that portion which is adjacent to the feasible 
region, R, for all of the constraints. 

At any particular i3 e R* it is useful to define the local 
boundary hypersurface, B (u) , to the complete feasible region as 
the intersection of the active constraints at 0. Let H(u) 
denote the set of indices of the tight constraints at Q. 

Then, symbolically, 

B(Q) - O B ± (VIII-21) 

icK(u) 

Clearly B(&) is an m-R dimensional nonlinear manifold in the 

— a 

m-dimensional independent variable space. 

An m-k^ dimensional linear manifold C(4) approximating 

B(d) is the intersection of the active linearised constraints at 
u; that is. 
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L 



t 


V ' 


T 


:(u) « n ^(u) 

icK(Q) 


|u:S(Q) (u - Q) + • o| 


(VIII-22: 

(VIII-23; 


Now let Q(G) denote the linear space spanned by the gradients 
to the active constraints; that is. 


Q(u)“|‘i : 3 <* v •••»«„ for which u - ^ Oj ,J (VIII- 


24) 




1-1 


and let Q(u) denote the orthogonal complement to Q(Q); that 
is. 


R® - Q(d) 0 Q(Q). 


(VIII-25) 


It can be shown that Q(Q) is the unique linear space that can 
be translated to obtain the linear manifold C(u). 


Furthermore there exist unique orthogonal projection oper- 
ators P(Q) and ^(u) that resolve any vector in the independent- 
variable space into its corresponding components in Q(u) arid 


\ jjj 

Q(u), respectively; that is, for any u c R 


u * p(u)u + P(fl)u, 


(VIII-26) 


where 


P(u)u t Q(u) and ^(u)u c Q(u). 
In particular. 


(VIII-27) 


P - S'(SS ') -1 s 


(VIII-28) 


and 


I - P, 


(VII1-29) 
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An additional concept is the idea of problem scaling. The 
purpose of problem scaling is to Increase the efficiency of the 
targeting/optimisation algorithms by transforming the original 
problem into an equivalent problem that is numerically easier to 
solve . 


To numerically scale a problem, two general types of scaling 
are required: (1) independent-variable scaling, and (2) dependent- 

variable scaling. Independent-variable scaling is accomplished 
by defining a positive diagonal scaling matrix, W^, such that, 

the weighted independent variables are given by 


% 

u * 



(VI 11-30) 


Simularly, dependent-variable weighting is accomplished by 
defining an optimization-variable scale factor, W p > an d * 

positive, diagonal, dependent -variable scaling matrix, W , such 
that the weighted optimization variable is e 


P x * W p F(u) 


(VIII-31) 


and the weighted dependent Variables are given by 


£(u) 



yielding a weighted error fun jt ion 


p 2 = (>)• 


(VIII-32) 


(VIII-33) 


The program contains several options for computing the in- 
dependent-variable weighting matrix. However, the option most 
often used is the percentage scaling matrix 


The dependent-variable weighting matrix is always computed 
as the reciprocal of the constraint tolerances, and is given by 

| W 1 - — , (VIII-35) 

L Mu e t 

th 

where is the tolerance for the i constraint. The optimiza- 
tion scale factor is merely Input so that is approximately 

equal to one. 

For simplicity, the following discussion of the algorithm 
assumes an appropriately scalec problem. However, the scaled 
equations can be obtained by making the following simple sub- 
stitutions: 

a 

u replaced by u 

F replaced by 

c replaced by £ 

b replaced by 

S replaced by |wj [S] ^ 

7F replaced by Wp W ^ VF. 

The final key concept employed by PGA is the idea of a direc- 
tion of search. Heuristically , the direction of search is nothing 
more than a particular line in the independent-variable space 
along which the constraint error is reduced, or along which the 
cost-function is decreased. In a more precise sense, the direc- 
tion of search at Q is a half-ray emanating from u. Thus, for 
any positive scalar, y, the equation 

u - fl + yi (VIII-36) 

sets the limits of this half-ray and represents '’movement 1 * in the 
direction a from u. This is illustrated in figure VIII-1. 
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Figure VIII-1.- Direction of Search in the 
Independent -Variable Space 

If § is a unit vector, then y represents the actual dis- 
tance "moved" in the direction £. This concept of direct ion-of- 
search is particularly important since it enables the m-dimen- 
sional nonlinear programming problem to be replaced by a sequence 
(hopefully finite) of one-dimensional minimizations. What rema 
to be explained then is: (1) how to select 

search; and (2) how to determine the step size in that directs . 
All "direct" optimization methods employ this concept and, hence, 
differ only in their answers to the two preceding questions. The 
technique by which s^ and are selected by PGA will be de- 

scribed in subsequent sections. 


Direction of Search 

The projected gradient method uses two basic search direc- 
tions. For this discussion, they will be termed the constraint 
and optimization directions, respectively. PGA proceeds by tak- 
ing successive steps in one or the other of these two directions. 
The computation of each of these search directions is described 
below at a particular point u in the independent-variable space 
where fi of the constraints are active. 
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Constraint direction .- The constraint direction depends 
critically on the number of active constraints. Three cases are 
distinguished below: 

1) Case 1.- If fi < m, then that unique control correc- 
a 

tion Afk is sought, which solves the linearized con- 
straint equation 

S(d) Aii + * 0^ (VIII-37) 

and minimizes the length of Au. The solutions to the 
preceding vector equations define the m-n^ dimensional 

linear manifold C(G), which approximates the local 
boundary at 0 as described in detail in the preceding 
section. The - desired minimum norm correction, AO, is 
then the vector of minimum length in the indepdendent- 
variable space from Q to the linear manifold C(0). 
Analytically, it is given as 

Ad - — S'lSST 1 *^). (VIII -38) 

This correction is illustrated in figure VIII-2. 

The direction of search then is simply taken to be this 
minimum-norm correction to the locally active linearized 
constraints; that is, 

s C (fi) - Ad. (VI IT -39) 



figure VIII-2.- Illustration of Minimum-Norm Constraint, 
Direction for ■ 2 < m ■ 3 
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. — ^ C(0) , intersection of 

First linearised linearis'd constraints 

constraint 

Figure VIII-3.- Illustration of Newton-Raphson Constraint, Direction 
for n # ■ a ■ 3 

The direction of search is taken to be this unique .or- 
rection vector satisfying the linearised constraints; 
that is. 


s c <d) ■ Ad. 


(VI1I-41) 
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is minimized with respect to Au. Gauss demonstrated 
that the formula for this "least squares" correction is 


AG - -(S'S)- 1 S'e(G) . 


(VIII-43) 


Figure VIII-4 illustrates the least-squares correction pi~~ 
torially. As in the preceding two cases, the search 
direction is then taken to be this optimal correction; 
that is. 


s C (G) - AG. 


(V1II-44) 


Third linearized 
constraint 


Fourth linearized 
constraint 


Second linearized 
constraint 


Second linearized 
constraint 



AG, least-squares 
correction 


First linearized 
constraint 


Figure VIII-4.- Illustration of Least-Squares Constraint, 
Direct for n # - 4 > m » 3 
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Optimization direction .- When the number of active constraints 
is less than the number of independent variables, it is then pos- 
sible to reduce the nonminimal cost-function. Obviously the 
steepest descent direction, -VF(fl), would be the best local 
search direction for reducing the cost function. Such a direc- 
tion, however, would generally produce unacceptable constraint 
violations. To avoid this difficulty PGA orthogonally projects 
the unconstrained negative gradient, - VF (u ) , into a direction 
parallel to the local linearized constraint boundary C(fl). By 
searching in the direction of this negative-projected gradient 
the algorithm can guarantee that there is no further constraint 
violation than that of u for the case of linear constraints. 

To calculate this direction, it is only necessary to apply to 
the unconstrained negative gradient the projection operator P(u), 
which maps any vector in the independent-variable space into its 
component in Q(u), the unique linear space that can be trans- 
lated into coincidence with the linear manifold C(u). Thus, 

s°(0) - -P(ti) VF(d) 

- -[I - F(fl)] VF(fl) 

- -[I -S' (SS')" 1 S(d)] VP(fi) 

The direction of search for the accelerated projected gradient 
method is 


(VIII-46) 


(VIII-47) 


s°(0) - -H P 7F(G) 
— n ^ n — — 

where 

H o- 1 


(VI 11-45) 


and 


H „ ■ H „-l + A n + V " h,t * " ' 2 

*n ‘ [ 4 ^ V 

»„ - -[«„-! lAVlJ/ft-A’ 

'in * Wl 1 


1 (VIII-48) 
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Figure VIII-5 illustrates the direction of the negati\ e-pro jected 
gradient for the case of a single active constraint. 


First Inset lv* 
const raint 


Linearized approximation 
to sole active constraint 

at i 

Unconstrained negative 
gradient to cost 
funct ion at i 


Negative protected 
gradient at u 



Sole act! '• 
constraint 
(paraboloid 
of revolution) 


Second Inactive 
constraint 


Figure VIH-5.- erection of Negative-Projected Gradient for n ft - 1 

and m - 3 (Feasible region is that region inside 
paraboloid, above lower plane, and below upper plane; 
cost— function is vertical height) 

If there are no equality constraints, and if all the inequality 
constraints are inactive, then S is the zero matrix and the 
direction of search becomes the standard deflected gradient 
direction 


s°(u) - -H n VF(G). 


(VI I 1-49) 


Similarly, if the single-penalty-function method* are used, 
then the directions of search that minimize 


P 2 . p + W e'e 


(VII 1-50) 


are: 

1) Steepest-descent method 

s°(d) - -VP 2 (fl); 
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2) Conjugate gradient method (steepest— descent starter) 

s , , where n > 2, 
— n-l 


s = -VP, u + 
— n — 2 — n 


' CM 
04, 
>1 

u 

lH 

i)— 

M 

1? 2< 

u 

l-n- 

■l) 

2M 

vi) 


3) Davidon's method (steepest-descent starter) 

s° * -H VPju ’ where n _> 2 
— n n — 2\— nj 


and 

H = H 
n 


1 i + A + B , 
n-1 n n 


where A and B have the same definitions as in the 
n n 

accelerated projected gradient mode. 


Step-Size Calculation 

At any particular point u in the independent— variable 
space, the PGA algorithm proceeds by reducing the multidimen- 
sional problem to a one— dimensional search along the constraint 
direction to minimize the sum of the squares of the constraint 
violations, or along the optimization direction to minimize the 
estimated net cost-function, in either case, once the initial 
point u and the direction ofsearch si are specified, the prob- 
lem reduces to the numerical minimization of a function of a single 
variable — namely , the step size. PGA performs this numerical 
minimization via polynominal interpolation, based on function 
values along ttie search ray and the function’s value and slope at 
die starting point. Consider then, in detail, the calculation of 
this latter pair of quantities for the respective functions asso- 
ciated with the constraint and optimization directions. 

Constraint direc tion .- The function to be minimized along 

the constraint direction, s_ C , is the sum of the squares of the 
const ra i nt v io ia t Lons; namely 

i, c <,> - (vi 1 1 -si) 


vm-16 


Clearly 

h (0) - |e(u) | 2 . 
c 

Differentiation via the chain rule yields 

h '(0) = 2e'(a)S(a)s C . 

C 


(VIII-52) 
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Recall that the search direction s was obtained as an in- 
dependent-variable correction either satisfying «11 the linearized 
constraint equations if n fl 1 m, or minimizing their violation 

n < n . Thus, if the constraints are reasonably linear, a good 
a 

initial estimate for the > minimizing h £ is one. 

Optimization direction .- The function to be minimized along 

the optimization direction, s°, is the estimated net cost- 
function which is defined as 

h (y) - F (u + ,s°) - F(u) + l'nu) [-S'(SS')- 1 e(u + ,s°)] . (VIII- 

O ' 7 t 

v . J v ' 
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change in cost- 
function produced 
by step of length 

, along £° 


linearized approximation to 
change in cost-function re- 
quired to perform minimura- 
norm correction back to the 
feasible region 


Clearly 


h (0) - -V/F(0) S *(SS ') -1 (u)e(u) . 
o 


By expanding h in a Taylor series in , about 


(VII 1-55) 


0 . 


and by making use of the fact that Ps° - 0 since s lies in 
Q(u), it can be shown that 

h '(0) - i°- (VIII-56) 


These properties of h Q are illustrated in figure 25. 
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Figure VIII-6.- Properties of Estimated Net Cost Function 

Both the constraint and optimization directions are based on 
a sensitivity matrix that depends critically on which constraints 
are active. Hence, for searches in either direction, it is im 
portant to limit the step size so that the set of active ® 

does not grow. Such a limit can be obtained based on linear ap 
proximation and suffices to deal with inactive constraints becom- 
ing active. 

The reverse situation— of active constraints becoming in- 
active-poses no difficulty. To see this, note that because of 
our treatment of the active constraints as linear manifolds, a 
first-order approximation of the distance to a particular active 
constraint boundary would not change along the optimisation direc- 
tion. Furthermore, along the constraint direction any change in 
the status of an active constraint will be appropriately treated 
by minimizing h with respect to the step length. 
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denote the set of active constraint indices at 


Let K(u) 
fl, and let 

, (VIII-57) 


where s(u) is the search direction at vector u. Then assign 
to each k in K the number 


-c k (0)/r k if r k 


< 0 


X(k) 


, R if 


r, >0 
k — 


I 

I 


(VIII-58) 


where R is a very large real number. Then A(k) is a linear 
approximation to the distance along the search ray from u to 

the boundary, B k> of the k Ch constraint. Hence a resonable 

upper bound for the step length is 

X - min [X(k)]. (VIII-59) 

kcK 


One-Dimensional Minimization 

Monovariant minimization in PGA is performed exclusively by 
polynominal interpolation. First the actual function, f, to be 
minimized is fitted with one or more quadratic or cubic poly- 
nomials until a sufficiently accurate curve fit, p, is ob- 
tained; that is. 


n 

p ( Y ) -^^a 1 Y 1 ^f(y) for all y of interest. 
i-0 


(VII1-60) 


Then the independent variable value, y", that minimizes f is 

approximated by the value, y”, which minimizes p. Clearly, 

y® can be determined analytically if n ^ 3. 

P 
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The minimi zat ion routine makes ingenious use of all the in 
formation it accumulates about f to obtain a good curve fit. 
First, f is fitted with a quadratic polynominal, Pi, based 
on: 

, 1 ) f( 0 ) 

i 2 ) f'( 0 ) 

j 3 ) w bere y* > 0 is an initial estimate of the 

y value that minimises f. 

| The coefficients of this quadratic polynominal are then calcu- 

i la ted from the formulas: 

a 0 - f( 0 ) 

ax - f'( 0 ) 

*2 - {*(£) - -]/*: + *i/v 

j The value of the Independent variable that minimizes this poly- 

nominal is 


| (VII1-61) 


y 7 - -aj/2 a 2 . (VIII-62) 

If y* and Y* do not differ significantly, y* is taken 
to be y® and the minimisation procedure is considered complete. 
Similarly, if Pi(y?) is not significantly different from 

f (y"), then y" is taken to be equal to y® and the process 
is terminated. Otherwise f is fitted with a cubic polynominal, 
p 2 , based on 

1^ f(0) 

2) f'(0) 

3) f (yo) and yo > 0 

4) f (y?). 


VI 11-20 


1 


If f is fitted using p 2 , then coefficients ere calculated 
from the following formulas: 

*\ 

a 0 * f(0) 


a 2 - f'(0) 

X - max (yS.yT) 
a - min (yo»Y?)^ 

a 2 - [Xaj a + a 0 (1 + a) + (a 2 f(*) - f(aX))/(l - o)]/(X 3 a 2 ) 
a 3 - [ (f (oX) a 3 f(X))/(l - a)-Aa(l + o)ai-(l + a - a 2 )a 0 ]/(X 2 o 2 )J 


V (VIII-63) 


The value of the independent variable, X 2 ,- that minimizes this 
cubic polynomial is 


Y2 " ( _ a 2 + v 4 2 ~ / 3a 3* 


(VIII -64) 


If y” and Y? do not differ significantly, y“ is taken 
to be y“ and th * minimization is stopped. Similarly, if P 2 (y 2 ) 
is not significantly different from f(Y 2 .), then y 1* taken 
to be equal to Y™ and the procedure is terminated. 


If none of these stopping conditions is met, a third quad- 
ratic curve-fit is attempted. The accumulated set of sample 

points on ' f , namely [0,f(0)], £yo» f ( Y o)]» |_Y?»^(y?)J» « nd 

Ty”» f (f“) 1 » !■ arranged in the order of their ascenditfg abscissa 

Values . Then the first point whose ordinate value is less than 
that of the following point is selected. 

To simplify the notation in the following pages, relable this 
point as Iy 2 , f(Y2>l. tha preceding point as Iyi. f(Yi)l. •«<* 
the following point as [Y3» f (y 3) ] * 
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Another quadratic polynomial, P 3 , is then fitted to 


1) f(n) 

2) f(Y 2 > 

3 ) 5 ( 0 . 


The formulas for these quadratic coefficients are as follows: 




Vj 


C ij " Y i + Y j 




y ±' y i 


ao 


b 2 3 bis bj 2 

d 12 d 13 1 d 21 d 23 d 31 d 32 

C 2 3 C 13 C 12 

ai " dTIdTI f(Yl) ' ^17 f(Y2) " dTidal f(Y3> 

a 2 ■ T^J— f CYl) + a K— f<Y 2 ) + 7~V~ f(Y 3>- 

2 d 12 d 13 d 21 d 23 d 31 d 32 
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The value of the independent variable that minimizes this quad- 
ratic is 


- -aj/2a 2 . (VI I 1-66) 

If y® and y" do not differ significantly, Y ffl is taken 
to be y™ and the search is discontinued. On the other hand, if 
p 3 (y®) is not significantly different from f(y?)» then y™ 
is taken to be (y*) end the process is terminated. 

If neither of these stopping conditions is met, then a cubic 
polynomial is fitted to 
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1) f(Yl). Yl 

2) £(y 2 ) » Y 2 

3) f(Y3>» Y 3 

4) f(Y4>» Y4 * Y?. 

The formulas for these coefficients are as follows! 

D 1 = (Y 2 ~ Yl)(Y3 ~ Y 1 ) <Y 4 “ Y.l) 

D 2 - (Yl ~ Y2>(Y3 ~ Y 2 ) (y 4 “ Y 2 ) 

D 3 “ (Yl ~ Y 3) (Y2 ~ Y3MY4 “ Y3) 

D4 - (y 1 - Y 4)(Y3 _ Y4) (Y3 ~ Y4) 

Y 2 Y 3 Y 4 Y 1 Y 3 Y 4 YiY 2 Y4 # % YlY 2 Y3 ^ ^ 

a 0 - — — f(Yi) + — £ (Y 2 ) + lx." f <Y3) + rT f (Y4> 


Di D 2 

Y2Y3 + Y2Y4 + Y3Y4 

Di 


d 3 


d 4 


(YlYi + Yl *4 + Y4Y3> 

f(Yi) + f (Y2 


D 2 


- 


I : 

'A « 


<YlY2 + Y1Y4 + Y 2 Y 4 > (YlY 2 + Y 1 Y 3 + Y2Y3) 

+ f (y 3) + ~ f (y 4) 


D3 

(Y 2 + Y 3 + Y41 


(Yl + Y3 + Y4) 


»2 


f(Y i> + d 2 

(Yl + Y2 + Y4) 


f(Y3) 


D4 
f(Y 2 ) 

(y 1 + Y 2 + Y4) 


f (Y4) 


a 3 - 


- — f (Yl) - ^ f (Y 2 ) - ^ £ (y 3) - ^ f (* 0 - 
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The value of the indpendent variable minimizing this fourth cubic 
polynomial is 


y® ■ (-a 2 +J*2 - 3* 3*1)/ 1*3' 


(VI1I-68) 
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If y“ and y“ do not differ significantly, Y® is taken 
to be y® an d the minimization is stepped. Similarly, if 
p 4 (y®) is not significantly different from f (y“) » then y® 
is taken to be equal to y® and the procedure is terminated. 

If none of these stopping conditions is met, the accumulated 
set of sample points is searched for the point with the minimum 
ordinate value. The abscissa value of this poiut is taken to be 

Y®, and the minimization is considered complete. 


Algorithm Macrologic 

After being initialized the projected gradient algorithm 
proceeds as a sequence of iterations, each consisting of an op- 
timization step followed by a constraint-correction step (see 
fig. VIII-7) . The very first 3 tep from the user s initial indepena 
ent-variable estimate is however, one of constraint correction. 
Furthermore, the optimization step is also omitted on any itera- 
tion for which the constraint-violation function, h £ , was not 

reduced by the constraint correction step of the preceding itera- 
tion. 


The optimization search direction that emanates for is 

based on the sensitivity matrix, S (^n) ; is# 

<vm - w 

as discussed previously. Hence, s£ lies in the subspace 

o 

The value of the independent-variable vector, u^ after 
the optimization is 


n 


* U + Y, 


where y is the optimum step size, 
o 


(VIII-70) 
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Figure VIII-7.- Macro logic of Projected 
Gradient Algor it ha 
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APPENDIX A 

DECOMPOSITION BY PARTITIONING INTO FULL-RANK SUBPROBLEMS 


Consider a trajectory consisting of s mission segments each of 
which may consist of one or more phases. Suppose that each seg- 
ment has its own physical control vector , u , containing m^ com- 
ponents and its own constraint vector c k having n fe components. 

Let v k be the target value of the constraint vector for segment 
k. Finally, suppose there are n constraints C^, which are best 

associated with the mission as a whole rather than any particular 
segment. Let denote the target values of these constraints. 

The problem is then 


minimize: 

s 




k 1 


(v •■)•£' '“(«■■) 

\k«l ' k-1 U-l / 


subject to: 


k k 
v i ” c i 


V i" C i 


te*') 

te-‘) 


for k“l, . . . ,s 
for i-l, • . . ,n. 


for i*l,...,n. 


When numerous segments are present with varying degrees °* * n ~ 
fluence on the overall objective, F, solution of problem [1J in 
a single piece by means of any existing equality-constraint 
minimization procedure becomes impractical if not impossible. 

By solving coordinated sets of subproblems representing the in- 
dividual segments, the decomposition procedure is able to solve 
problem 11] routinely. The decomposition technique thus both 
imitates and extends the intuitive approach of the experienced 
trajectory designer. 


[ 1 ] 


The notation -U- denotes the ordered union of the indexed quantity 

A-l 


k«l 

immediately to its right 


The process of decomposing a problem by partitioning it into full- 
rank subproblems is based upon two fundamental ploys • The first 
is the use in each segment of certain key control variables to 
satisfy the constraints of that segment. Thus segment k is made 
into a full-rank subproblem by designating n^ of the physical con- 
trol vector components as subproblem variables and using them 
to satisfy the constraints of that segment. The remaining 

control vector components of segment k are grouped together with 
similar variables from the other segments. This collection of 
controls, together with the overall mission objective, F, and 
constraints, £, are made into a master problem of minimization 
subject to equality constraints. This partitioning of the orig- 
inal problem lends itself well to computation. The subproblems, 
which must be re-solved for each new choice of master problem 
controls, are solved by the highly efficient Newton-Raphson pro- 
cedure. The master-problem, which serves to coordinate the sub- 
problem solutions, uses the equally efficient but more time-con- 
suming accelerated projected-gradient algorithm. 

The second fundamental ploy is use of constraint target values 
of the various subproblems as ma3ter-problem independent variables. 
To obtain an optimum composite trajectory from a set of mission 
segments, the mission analyst typically varies the segment aim 
points parametrically and chooses the endpoint combination that 
results in the lowest overall cost. Indeed, in trajectory 
analysis the decision "where to go” is usually more important 
than "how to get there." By using subproblem constraint targets 
as master-problem controls the decomposition procedure automates 
the analyst's successful design approach. 

The successful convergence of the decomposition procedure demands 
a reasonable partitioning of the original controls and constraints 
into the master-problem and subproblem categories. To be more 
precise , the subproblem controls and constraints must be so chosen 
that each subproblem will have a solution for any sat of subproblem 
constraint target values that might reasonably arise during the 
master-problem iteration process. Thus, for a given subproblem, 
the controls that have the most substantial effect on that subprob- 
lem’s constraint set should be chosen. Similarly if a particular 
constraint cannot be assigned to any subproblem whose controls 
can achieve its satisfaction, it should be designated a master- 
problem constraint. Finally, the number of master-problem con- 
strains should be held to a minimum. Indeed the master-problem 
should be kept as simple as possible because each of its iterations 
requires the re-solution of all the subproblems. 
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The decomposition procedure maintains simulation flexibility in 
obtaining master-problem control sensitivity information by using 
numerical differencing. Solution of the master-problem by descent 
requires constrained derivatives — quotients of dependent perturba- 
tions in the master-problem objectives and constraints by inde- 
pendent perturbations in the master-problem controls assuming tb&t 
the subproblem controls adjust uniquely to keep the subproblem 
targeted. These derivatives could be approximated by the numerical 
differencing of master-problem trajectories consisting of iterative- 
ly targeted subproblems. This approach, however, must be rejected 
because it is both susceptible to numerical error and demanding 
in computational effort. Instead formulas are used that relate 
the constrained derivatives to the partial derivatives of the 
master-problem objective and constraints with respect to all of 
the physical controls of the original problem. These partial 
derivatives are approximated conveniently by numercial differ- 
encing of the master-problem trajectories without subproblem 
targeting. Thus, the need for depriving variational equations 
for each simulated trajectory is eliminated for a reasonable 
computational price. Any trajectory that can be simulated, can 
be shaped with no additional analytical effort. 


To precisely define the procedure, considerable nomenclature 
must be established. Most of the user supplied parameters have 
already been defined. Two, however, remain. The first is p^, 

the number of subproblem constraint target values from segment 
k which are to be used as master-problem independent variables. 
The second is q^, the number of master-problem constraints aris- 
ing from segment k. 


Consider next the procedure's working variables. To simplify 
notation, the subproblem and master -problem controls are given 

distinct literal symbols and resequenced. Indeed, let y^ denote 

tne jth subpnblem control arising from segment k and be the 

jth master-problem control from that segment. Similarly, the 
segment constraint target values are assigned new symbols to 
distinguish those that are to be held fi^ed from those that are 
to be used as majter-problera Independent variables. Indeed, let 
v 

r. denote the ith constraint target value from segment k that 

1 k 

is to be used as a master-problem independent variable, and t^ , 

the 1th constraint target value to be held fixed for that segment. 
Finally, the mas ter -problem constraint vector is resequenced so 
that Its first qi components denoted by arise from the first 
subproblem, and the next qj components denoted by C 2 arise from 
the second subproblem, and so forth. In terms of this new nota- 
tion, the original problem (1] becomes 


/ 


minimize : 




y ( i 1 ' » * k ) 


, (), k 0 * k ) 


k-1 
subject to: 
s 


\u 

k*l 


U u ‘1 


k=l 


where : 


m » *1 


f OX 1 ) • • • f 8 


I k-1 


£ is a master-problem independent variable 

lc 

£ is a master-problem independent variable 


- y - 

k=l 


(r k U i k ) 


is the unique vector of subproblem independent 


variables for subproblem k, satisfying that subproblem's constraint 
set for the current set of master-problem independent variables. 


To solve problem [2], a procedure must first be established for 
solving the subproblems , that is, for determining the unique 
k k k 

X given the £ and £ • As noted above, the Newton-Raphson 
algorithm for solving full-rank systems of nonlinear equations 
is the technique selected. To start the iterative solution, the 

user must input a good estimate, (x k ) Q i for the physical-control 

vector of subproblem k, which approximately yield the constraint 
target values for that subproblem. The procedure, then, succes- 
sively refines this estimate using the Newton-Raphson recursion 
formula 





v 


where 

«k 3s - 

A • — for i»l, . .. ,s 
for k*l, . . . 
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After the first set of subproblems are solved, the converged sub- 
problem control values from one subproblem set are used as the 
starting estimates for the next subproblem set* lurther, the 

kk . . 

Jacobian matrix; A , is updated from one iteration to the next 

only if the old Jacobian does not reduce the constraint errors 

in norm by a user specified fraction, p. 

The equality constrained minimization that is the master -problem 
is carried out by the projected gradient algorithm already familiar 
to POST users. The only new technique involved is the computa- 
tion of the constrained derivatives in terms of the partial 
derivatives of the master-problem objective and constraints with 
respect to all of the physical controls of the original problem. 
First, the perturbation of the subproblem independent variables 
caused by perturbations in the master-problem independent vari- 
ables must be determined. Both master-problem physical con- 
trols and subproblem constraint target values must be consid- 
ered. Once these constrained derivatives of the subproblem 
independent variables are determined , th>y can be used to 
calculate the desired constrained derivatives of the master- 
problem objective and constraints with respect to all of the 
master-problem independent variables. 

The constrained derivatives of the subproblem controls are all 
derived from the basic subproblem equation 



2 

r 


2.-1 

£ 

k«l 


2k k 
A i 


B 


2k k 


k"l 


for 2*1 , . . . ,s 


where 


.2k 



for 2*1, . . . ,s. 
for k*l,. . . ,2 


The constrained derivatives with respect to the master-problem 
physical controls are given by the equation 
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o-k 



for f*l 
for k*l 


,8 
2 . 


[ 5 ] 


[ 6 ] 


» * • • * 


17 ] 


The constrained < erivatives with respect to the subproblem con- 
straint target values are computed from the two formulas 


ilz; 


- (* kk ) 


for k=l s 


[ 8 ] 


P k . 


where is the matrix consisting of the first p R columns of 


the identity matrix of order n, , 

k 




■.-I 


6 jl 


o=k 


6r 


for U=l,. . . ,s. 
for k«l, . . . ,e-l 


[9] 


fht constrained derivatives of the master-problem objective with 
respect to both the master-problem physical controls and the sub- 
problem constraint target values follow from the ,, chain-rule" for 
differentiation. They are computed as 


■jF 


6z 


■ =k 


b ,k +' 


, o 
•o 


o=k 


for k=l , . . . 


[ 10 ] 


and 


s i 


6F 

6 


Z. n — u 1. ^ r 


for k=l, . . . ,s. 


Ul] 


— f=k o=k 
where 


• k 3f 


f. 


3 X 


for 1=1,, 
for k"l,, 


»s 


112] 


and 

. ' k 3f 4 


d Z 


for ■ ■ 1 , , 
lor k*l,, 


[13] 
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Finally, the constrained derivatives of the waster-problem con- 
straints with respect to the mas ter -problem physical controls 
and the subproblem constraint target values follow from a stra:ight- 
f rward application of the "chain rule." They are related to the 
appropriate partial derivatives by the equations 




for 1 , * • * , s 
for k*l , . . . , i * 


[14] 


and 


6C* A 6*° 


6r 


o=k 


6r 


where 


tk 3 £ 

> -71 

3z 


and 


_ * 




for £-1 , . . . »s 
for k=l , . . . , i 9 


for . ,s 

for k+1 , . • . , l * 


for ^*“1. } • « • f s 
f or k=l # * 


[15] 


[16] 


[17] 
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